bookboon.com 


Boundary Element Methods for 
Engineers: Part II 

Plane Elastic Problems 
Roger Fenner 



Download free books at 

bookboon.com 

















































































































































































Roger Fenner 


Boundary Element Methods 
for Engineers 

Part II: Plane Elastic Problems 


2 


Download free eBooks at bookboon.com 



Boundary Element Methods for Engineers: Part II: Plane Elastic Problems 
1 st edition 

© 2014 Roger Fenner & bookboon.com 
ISBN 978-87-403-0733-7 


3 


Download free eBooks at bookboon.com 



Boundary Element Methods for Engineers: 
Part II: Plane Elastic Problems 


Contents 


Contents 



Preface 

Part 


Notation 

Part 


Some Program Variable Names 

Part 

1 

Introduction 

Part 

1.1 

Continuum Mechanics Problems 

Part 

1.2 

Some Practical Engineering Problems 

Part 

1.3 

Methods for Solving Harmonic and Biharmonic Equations 

Part 

2 

Boundary Element Analysis of Potential Problems 

Part 

2.1 

Fundamental Solution 

Part 

2.2 

Boundary Integral Equation 

Part 

2.3 

Discretisation of the Boundary Integral Equation 

Part 

2.4 

Constant Boundary Elements 

Part 



We do not reinvent 
the wheel we reinvent 


Fascinating lighting offers an infinite spectrum of 
possibilities: Innovative technologies and new 
markets provide both opportunities and challenges. 
An environment in which your expertise is in high 
demand. Enjoy the supportive working atmosphere 
within our global group and benefit from international 
career paths. Implement sustainable ideas in close 
cooperation with other specialists and contribute to 
influencing our future. Come and join us in reinventing 
light every day. 


OSRAM 

SYLVAN!A 


Light is OSRAM 


www.sylvania.com 


4 

Download free eBooks at bookboon.com 










Boundary Element Methods for Engineers: 
Part II: Plane Elastic Problems 


Contents 


2.5 

2.6 

2.7 

2.8 

3 

3.1 

3.2 

3.3 

4 

4.1 

4.2 

4.3 

4.4 

4.5 


Quadratic Boundary Elements Part I 

Scaling Part I 

Solving the Linear Equations Part I 

Solving Poissons Equation Part I 

Constant Boundary Element Program for Potential Problems Part I 

Program BEM2PC Part I 

Some Test Problems for BEM2PC Part I 

An Example: Downstream Viscous Flow in a Rectangular Channel Part I 

Quadratic Boundary Element Program for Potential Problems Part I 

Program BEM2PQ Part I 

Some Test Problems for BEM2PQ Part I 

An Example: Downstream Viscous Flow in a Rectangular Channel Part I 

An Example: Heat Conduction in a Domain of Complex Shape Part I 

Discussion Part I 

Appendix A: Gaussian Quadrature Part I 

Appendix B: Gaussian Elimination Part I 



Deloitte 




Discover the truth at www.deloitte.ca/careers 


© Deloitte & Touche LLP and affiliated entities. 


5 

Download free eBooks at bookboon.com 




Boundary Element Methods for Engineers: 
Part II: Plane Elastic Problems 


Contents 


Appendix C: Matlab Version of Constant Boundary Element 

Program for Potential Problems Part I 

Appendix D: Matlab Version of Quadratic Boundary Element 

Program for Potential Problems Part I 

Solutions to Problems - Part I Part I 

Preface 9 

Notation 10 

Some Program Variable Names 13 

5 Boundary Element Analysis of Plane Elastic Problems 20 

5.1 Fundamental Solution 21 

5.2 Boundary Integral Equations 32 

5.3 Discretisation of the Boundary Integral Equations 35 

5.4 Boundary Conditions and Surface Stresses 37 


SIMPLY CLEVER 


SKODA 



We will turn your CV into 
an opportunity of a lifetime 


Do you like cars? Would you like to be a part of a successful brand? 
We will appreciate and reward both your enthusiasm and talent. 
Send us your CV. You will be surprised where it can take you. 


Send us your CV on 

www.employerforlife.com 


6 

Download free eBooks at bookboon.com 












Boundary Element Methods for Engineers: 
Part II: Plane Elastic Problems 


Contents 


5.5 Quadratic Boundary Elements 41 

5.6 Scaling 51 

5.7 Solving the Linear Equations 51 

5.8 Body Forces 51 

6 Quadratic Boundary Element Program for Plane Elastic Problems 53 

6.1 Program BEM2EQ 54 

6.2 Some Test Problems for BEM2EQ 100 

6.3 An Example: Confined Compression of a Rubber Block 116 

6.4 An Example: Stress Concentration at a Hole in a Flat Plate 118 

6.5 Discussion 122 

7 Further Applications 128 

7.1 Axi-symmetric Problems 128 

7.2 Higher-Order Boundary Elements 128 

7.3 Three-dimensional Problems 128 

7.4 Non-Linear Problems 130 

7.5 Comparison with Other Methods 130 




MAERSK 


I joined MITAS because 
I wanted real responsibility 


The Graduate Programme 
for Engineers and Geoscientists 

www.discovermitas.com 


Real work 
International opportunities 
Three work placements 


0 


a 


I was a construction 
supervisor in 
the North Sea 
advising and 
helping foremen 
solve problems 



Download free eBooks at bookboon.com 









Boundary Element Methods for Engineers: 
Part II: Plane Elastic Problems 


Contents 


Appendix A: Gaussian Quadrature 132 

Appendix B: Gaussian Elimination 135 


Appendix E: Matlab Version of Quadratic Boundary Element Program 
for Plane Elastic Problems 


141 


Solutions to Problems - Part II 


188 




Because achieving your dreams is your greatest challenge. IE Business School’s Master in Management taught in English, 
Spanish or bilingually, trains young high performance professionals at the beginning of their career through an innovative 
and stimulating program that will help them reach their full potential. 

Choose your area of specialization. 

Customize your master through the different options offered. 

Global Immersion Weeks in locations such as London, Silicon Valley or Shanghai. 

Because you change , we change with you . 


www.ie.edu/master-management mim.admissions@ie.edu f # In YnTube ii 


Master in Management 


8 

Download free eBooks at bookboon.com 












Boundary Element Methods for Engineers: 
Part II: Plane Elastic Problems 


Preface 


Preface 


A few decades ago, the advent of high-speed electronic digital computers gave tremendous impetus to all 
numerical methods for solving engineering problems, and made it possible to solve with good accuracy 
many problems which previously could only be solved approximately. Finite difference methods, applied 
manually, gave way to finite element methods, which are still one of the most versatile and widely used, 
particularly in structural and solid mechanics. In thermofluids, methods of the finite volume type tend 
to be preferred. 

Slower to develop have been boundary element methods, based on boundary integral equations. Initial 
development was largely in the hands of mathematicians, as the underlying mathematics are relatively 
sophisticated. It was engineers, however, who turned boundary element methods into practically useful 
and powerful techniques. 

The purpose of this book is to serve as a deliberately simple introduction to boundary element methods 
applicable to a wide range of engineering problems. The mathematics are kept as simple as reasonably 
possible. Computer programs form an integral part of the boundary element approach and they are 
treated as such in the text. Several programs suitable for use on desktops or laptops are presented and 
described in detail and their uses are illustrated with the aid of a number of practical examples. Problems, 
with solutions, are provided at the ends of the chapters, for readers to solve for themselves. 

The programming language used in the main text is Fortran. Although it is somewhat unfashionable these 
days for general programming purposes, Fortran is still very widely used in engineering computation. 
Matlab versions of the programs are also provided in Appendices. Full listings of all the programs, both 
Fortran and Matlab, are available for download here . 

A prior knowledge of either Fortran or Matlab is desirable. The level of continuum mechanics, numerical 
analysis, matrix algebra, vector analysis and other mathematics employed is that normally taught in 
undergraduate engineering courses. The book is therefore suitable for engineering undergraduates and 
other students at an equivalent level. Postgraduates and practising engineers may also find it useful if 
they are comparatively new to boundary element methods. 

The book is presented in two Parts. Part I started with a brief review of the problems encountered in 
engineering, showing that they of two broad types. It then described boundary element treatments of 
problems of the potential type, using both constant and quadratic boundary elements. This Part II is 
concerned with elastic stress analysis problems of the plane strain and plane stress types. 

Imperial College London Professor Roger Fenner 
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Notation 


Notation 


The mathematical symbols commonly used in the main text are defined in the following list. In some 
cases particular symbols have more than one meaning in different parts of the book, although this should 
not cause any serious ambiguity. 


A 

[A] 

Ay 

d 1 ; 0*2, &3 

[B] 

B H 

[b] 

c 

c c c c 

u xx > u xy > u yx > u yy 

D 

E 


e 

A 

A 

G 

9 

9 

H 

H 

h 

i 

i 

J 

j 

j 

K 

k 

k 

k 

L 

M 

m 

N 

N c 

n 

n 


area of a solution domain 
square matrix 
coefficient of matrix [^4] 

constants in general boundary condition Equation 1.83 

square matrix 

coefficient of matrix [B] 

column vector of known coefficients 

torsional couple 

specific heat 

free term constants 

flexural rigidity of a flat plate 

Young's modulus 

strain 

drag and pressure flow shape factors for downstream flow 

function of position in Poisson’s equation 

function of position in biharmonic equation 

shear modulus 

acceleration due to gravity 

heat generated per unit volume 

height of a channel or solution domain in general 

total rate of heat conduction 

surface heat transfer coefficient 

nodal point number 

unit vector in the x co-ordinate direction 

Jacobian of transformation (from global to intrinsic co-ordinates) 
nodal point number 

unit vector in the y co-ordinate direction 
ratio of outer to inner radius for a cylinder 
thermal conductivity 
nodal point number 

unit vector in the z co-ordinate direction 

maximum dimension of the solution domain 

total number of boundary elements in a mesh 

element number 

total number of nodes in a mesh 

shape function 

direction of outward normal to the boundary of a solution domain 
outward normal vector to boundary 
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Tl X , Tty 

n 

Ti X , Tty 

P 

Pz 

v 

V 

Q 

Q 

R 

R 

r 

r 

r 

f 

S 

S 

s 

s 

9 

S 

T 

T 

T m 

t 

t 

U 

u 

u r 

u e 

u 

V z 

V 

V 

W 

w 

X,Y 
X, Y,Z 

x, y,z 

M 


components of n in the x and y directions 
unit outward normal vector to boundary 
components of n in the x and y directions 
source/force point on a boundary 
pressure gradient in the zco-ordinate direction 
pressure 

source/force point 
volumetric flow rate 
field point on a boundary 
field point 

a function of intrinsic co-ordinate <f 
radial co-ordinate 

distance between a source/force point and a field point 

vector distance between points 

unitvector in the direction between points 

boundary of a domain 

distance along a boundary 

vector tangent to a boundary 

ratio between lengths of successive elements on a boundary segment 

part of a boundary forming element m 

distance along a solution domain boundary 

temperature 

traction kernel function 

remote temperature of surroundings in thermal convection 

time 

traction 

displacement kernel function 
displacement or velocity in thexdirection 
displacement in the radial direction 
displacement in thexdirection 
mean value of u 

velocity component of a boundary in the z co-ordinate direction 
displacement or velocity in theydirection 
mean value of v 

width of a channel or solution domain in general 
displacement or velocity in thezdirection 
global Cartesian co-ordinates 

components of body forces per unit volume in the Cartesian co-ordinate 
directions 

Cartesian co-ordinates 

column vector of unknown quantities 


a coefficient of linear thermal expansion 

a coefficient in mixed boundary condition, Equation 2.32 

/? coefficient in mixed boundary condition, Equation 2.32 
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r 

AT 

£ 

V 

e 

e 

e 

K 

T 

V 

z 

7Tp 

TTq 

P 

a 

0 

0 

X 

0 

0 

V 
V 2 
V 4 


an angle 

change in temperature 
a small distance 

local intrinsic co-ordinate within an element 
angle of rotation per unit length of a bar in torsion 
angular co-ordinate 
an angle 

permeability of a porous medium 

viscosity 

Poisson's ratio 

local intrinsic co-ordinate within an element 

dimensionless pressure gradient 

dimensionless flow rate 

density 

stress 

velocity potential 

fundamental solution for potential 

stress function 

stream function 

potential 

grad operator 

harmonic (Laplacean) operator 
biharmonic operator 


Subscripts 

c 

e 

i,j,k 

L 

m 

n 

PI 

r 

s 

s 

T 

x, y,z 
6 

1,2 

1,2,3 


counter for nodes within an element 
von Mises equivalent (stress) 
nodal point numbers 
solution to Laplace’s equation 
element number 

direction of the outward normal to a boundary 
particular integral solution satisfying Poisson’s equation 
radial direction in polar co-ordinates 
direction along a boundary 
segment 

total solution (Laplace plus particular integral) 

Cartesian co-ordinate directions 

angular direction in polar co-ordinates 

inner and outer of two concentric circular boundaries 

nodes of a quadratic element 


Superscripts 

* effective value under plane stress conditions 

* modified quantity 
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Some Program Variable Names 

The Fortran computer program variable names widely used in the programs and main text are defined 
in alphabetical order in the following list. 


A Coefficients of matrix [A] 

All Matrix diagonal coefficient A..(potential problems) 

AIIXX, AIIXY, AIIYX, AIIYY 

Coefficients at the diagonal of matrix (elastic problems) 

AIJ Matrix coefficient A.. 

AK First kernel function contributing to the matrix [A] (potential problems) 

AKXX, AKXY, AKYX, AKYY 

First kernel functions contributing to the [A] matrix (elastic problems) 
Perpendicular distance from centre of curvature to mid point of a segment chord 
Square of ALPERP 

Element values of constants in mixed boundary conditions 

Nodal point values of constants in mixed boundary conditions (quadratic elements) 
Boundary segment values of constants in mixed boundary conditions 
Length of a segment chord measured between its end points 
Angular position of current end point on a curved boundary segment 
Angular position of first end point on a curved boundary segment 
Angle subtended at centre of curvature by a curved boundary segment 
Angular positions of end points on curved boundary segments 
Array storing element node contributions to the [A] matrix (potential problems) 

Array storing element node contributions to the [A] matrix (elastic problems) 

Array storing element node contributions to the [A] matrix (elastic problems) 

Coefficient of right hand side vector (matrix [B] times vector of knowns) 

Element values of constants in mixed boundary conditions 

Nodal point values of constants in mixed boundary conditions (quadratic elements) 
Boundary segment values of constants in mixed boundary conditions 


ALPERP 

ALPERP2 

ALPHA 

ALPHAN 

ALPHASEG 

ALSEG 

ANG 

ANGFIR 

ANGSEG 

ANGSTORE 

AROW 

AROWX 

AROWY 

BDPSI 

BETA 

BETAN 

BETASEG 

BIJ 

BK 


Matrix coefficient B.. 


Second kernel function contributing to the [B] matrix (potential problems) 

BKXX, BKXY, BKYX, BKYY 

Second kernel functions contributing to the [B] matrix (elastic problems) 

BK2 Non-singular part of second kernel function when P is the current element node 

(potential problems) 

BK2XX, BK2XY, BK2YX, BK2YY 

Non-singular parts of second kernel functions when P is the current element node 
(elastic problems) 
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BROW 

BROWX 

BROWY 

BTX, BTY 

CASE 

D 

DPSI 

DPSIPI 

DPSIPIM 

DPSISEG 

DPSISTORE 

DPSIT 

DRDN 

DUDZ 

DVDZ 

DZDE 

E 

EGL 

ELENGTH 

ESS 


Array storing element node contributions to the [B] matrix (potential problems) 

Array storing element node contributions to the [B] matrix (elastic problems) 

Array storing element node contributions to the [ B ] matrix (elastic problems) 

Coefficients of right hand side column vector (matrix [B] times the vector of knowns) 

Alphanumeric plane stress or strain problem type 

Perpendicular distance from P to the element containing node Q 

Nodal point values of the potential gradient solution to Laplaces equation 

Nodal point values of the particular integral potential gradient function 

Values of the particular integral potential gradient at the nodes of each element 

Values of potential gradient applied as boundary conditions to the boundary segments 

Temporary store for potential gradient 

Nodal point values of total potential gradient (Laplace plus particular integral) 

Rate of change of radius with distance along normal to boundary 
Rate of change of displacement u with £ along element 
Rate of change of displacement v with £ along element 
Jacobian of transformation from intrinsic co-ordinate £ to r\ 

Youngs modulus 

Values of the intrinsic co-ordinate at the Gauss points (logarithmic quadrature) 
Lengths of the elements 
Direct strain along boundary 
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ESTORE 

ETA 

Stored value of Youngs modulus 

Intrinsic co-ordinate r\ 

EVX 

x component of the vector along an element 

EVY 

y component of the vector along an element 

FI 

Constant functionin Poissons equation 

FLOWELEM 

FLOWIN 

FLOWOUT 

Potential flow across an element 

Total potential flow into the domain 

Total potential flow out of the domain 

FLOWSEG 

Flows of potential across the boundary segments 

FXELEM 

Force on an element in v direction 

FXSEG 

Total force on a boundary segment in v direction 

FYELEM 

Force on an element in y direction 

FYSEG 

Total force on a boundary segment in y direction 

HX 

Interval between points in the v direction used in domain integration 

HY 

Interval between points in the y direction used in domain integration 

I 

Node counter 

11,12,13 

IBC 

Numbers of the three nodes of a quadratic element 

Type number of boundary conditions applied to the (constant) elements 

IBCD 

Counter for segments subject to applied potential gradient boundary conditions 

IBCE 

IBCM 

IBCN 

IBCP 

Type number of boundary conditions applied to the (quadratic) elements 

Counter for segments subject to applied mixed boundary conditions 

Type number of boundary conditions applied to the nodes (of quadratic elements) 
Counter for segments subject to applied potential boundary conditions 

IBCPC 

Counter for point displacement constraints 

IBCS 

IBCU 

IBOUND 

IC 

IDIRPC 

Counter for segments subject to applied stress boundary conditions 

Counter for segments subject to applied displacement boundary conditions 

Counter for boundaries 

Case number for logarithmic Gaussian quadrature 

Direction numbers of point displacement constraints 

IEEND 

Counter for element end points 

IEP1 

Counter for first end point of an element 

IEP2 

Counter for second end point of an element 

IFIRST 

IFLAG 

Numbers of first nodes on the segments 

Flag for ill-conditioning of the [A] matrix 

IGAUSS 

IINT 

ILAST 

Counter for Gauss points 

Counter for internal points 

Numbers of last nodes on the segments 

IN 

Counter for nodes within an element 

IROW 

Number of row in the [A] matrix 
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ISEG 

ISEGBC 

ISEGELEM 

ISEGEND 

ISEGMAX 

ISEGMIN 

ISEND 

IT 

IX 

IXMAX 

IY 

IYMAX 

J 

JACOB 

JMAX 

M 

Ml 

M3 

MAXL 

MAXNB 

MAXNEL 

MAXNEQN 

MAXNNP 

MAXNPC 

MFIRST 

MLAST 

MMAX 

MMIN 

NBCD 

NBCM 

NBCP 

NBCPC 

NBCS 

NBCT 

NBCU 

NBOUND 

NEEND 

NEL 

NELB 


Segment counter 

Segment numbers for a particular type of boundary condition 

Segment numbers for elements 

Segment numbers for element end points 

Number of last segment on current boundary 

Number of first segment on current boundary 

Counter for boundary segment end points 

Number indicating type of Gaussian quadrature (normal or logarithmic) 
Counter for points in the % direction used in domain integration 
Maximum value of IX 

Counter for points in the y direction used in domain integration 
Maximum value of IY 
Node counter 

Jacobian of transformation from global to local intrinsic £ co-ordinate 
Maximum number of columns in the extended [A] matrix 
Element counter 

Numbers of the elements adjacent to the first node of each element 

Numbers of the elements adjacent to the third node of each element 

Maximum dimension of the solution domain 

Maximum number of boundaries allowed by the array dimensions 

Maximum number of elements 

Maximum number of equations 

Maximum number of nodal points allowed by the array dimensions 

Maximum number of point displacement constraints 

Numbers of the first elements on the segments 

Numbers of the last elements on the segments 

Number of last element on current boundary 

Number of first element on current boundary 

Number of segments subject to applied potential gradient boundary conditions 
Number of segments subject to applied mixed boundary conditions 
Number of segments subject to applied potential boundary conditions 
Number of point displacement constraints 

Number of segments subject to applied stress boundary conditions 

Total number of segments subject to applied boundary conditions 

Number of segments subject to applied displacement boundary conditions 

Number of boundaries 

Number of element end points 

Number of elements 

Numbers of elements on the boundaries 
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NELSEG 

NEP1 

NEP2 

NEQN 

NGAUSS 

NINT 

NNP 

NNPB 

NODE 

NODEPC 

NSEGB 

NSEGTOT 

NU 

NX 

NY 

PI 

PSI 

PSIBOT 

PSIIP 

PSIIPT 


Number of elements on current boundary segment 

Numbers of the first end points of the elements 

Numbers of the second end points of the elements 

Number of equations 

Number of Gauss points 

Number of internal points 

Number of nodal points 

Numbers of nodal points on each of the boundaries 
Numbers of the nodes of the elements 

Numbers of nodes subjected to point displacement constraints 
Numbers of boundary segments on each of the boundaries 
Total number of boundary segments 
Poissons ratio 

Number of internal points in the v direction used in domain integration 
Number of internal points in the y direction used in domain integration 
n 

Nodal point values of the potential solution to Laplaces equation 
Value of potential on the bottom edge of a rectangular domain 
Laplace equation potential at an internal point 
Total potential at an internal point (Laplace plus particular integral) 
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PSILEFT 

PSIPI 

PSIRIGHT 

PSISEG 

PSIT 

PSITOP 

PSIVAL 

R1 

R1X 

R1Y 

R2 

R2X 

R2Y 

RATSEG 

RFN 

RSEG 

RX 

RY 

SD 

SDL 

SF 

SFL 

SFN 

SIGE 

SIGNN 

SIGNNSEG 

SIGSN 

SIGSNSEG 

SIGSS 

STORE 

TITLE 

TMX 

TX 

TMY 

TY 

U 

UELEM 

UNGX 

UNGY 

UNMX 

Value of potential on the left hand edge of a rectangular domain 

Nodal point values of the particular integral potential function 

Value of potential on the right hand edge of a rectangular domain 

Values of potential applied as boundary conditions to the boundary segments 
Nodal point values of total potential (Laplace plus particular integral) 

Value of potential on the top edge of a rectangular domain 

Values of potential stored for domain integration 

Distance from point P to the first end of element containing node Q 
v component of radius vector from P to the first end of element containing Q 

y component of radius vector from P to the first end of element containing Q 

Distance from point P to the second end of element containing node Q 
v component of radius vector from P to the second end of element containing Q 

y component of radius vector from P to the second end of element containing Q 

Ratio between successive element lengths on current boundary segment 

Value of function R (<f) 

Radius of curvature of current boundary segment 

Component in v direction of unit radius vector from P to Q 

Component in y direction of unit radius vector from P to Q 

Shape function derivatives for quadratic elements (normal quadrature) 

Shape function derivative values for quadratic elements (logarithmic quadrature) 
Shape function values for quadratic elements (normal quadrature) 

Shape function values for quadratic elements (logarithmic quadrature) 

Shape function value at a Gauss point 

Nodal point values of von Mises equivalent stress 

Nodal point values of direct stress normal to boundary 

Boundary segment values of direct stress normal to boundary 

Nodal point values of shear stress along boundary 

Boundary segment values of shear stress along boundary 

Nodal point values of direct stress along boundary 

Stored values in the boundary condition application process 

Alphanumeric title for the problem (maximum 80 characters) 

Element nodal point values of traction in v direction 

Nodal point values of traction in v direction 

Element nodal point values of traction in y direction 

Nodal point values of traction in y direction 

Nodal point values of displacement in x direction 

Element values of displacement in x direction 
v component of the unit normal at a Gauss point 
y component of the unit normal at a Gauss point 
v components of the unit normals at the nodes of each element 
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UNMY 

UNX 

UNY 

USEG 

UV 

V 

VELEM 

VSEG 

WG 

WGL 

XC 

XCENT 

XEEND 

XFIRST 

XINT 

XLAST 

XMID 

XNODE 

XP 

XPOINT 

XQ 

XSEND 

XX 

YC 

YCENT 

YEEND 

YFIRST 

YINT 

YINTGL 

YLAST 

YMID 

YNODE 

YP 

YPOINT 

YQ 

YSEND 

YY 

ZETA 

ZG 


y components of the unit normals at the nodes of each element 
x components of the unit normals at the nodes 
y components of the unit normals at the nodes 
Boundary segment values of displacement in * direction 
Nodal point values of computed displacements (or tractions) 

Nodal point values of displacement in y direction 
Element values of displacement in y direction 
Boundary segment values of displacement in y direction 
Values of the Gaussian weighting factors (normal quadrature) 

Values of the Gaussian weighting factors (logarithmic quadrature) 
v co-ordinate of the origin for the particular integral function 
v co-ordinate of the centre of curvature of a curved boundary segment 
v co-ordinates of the element end points 
v co-ordinate of first end point of current boundary segment 
v co-ordinate of an internal point 

v co-ordinate of last end point of current boundary segment 

v co-ordinate of the mid point between the ends of a curved segment 

v co-ordinates of the nodes 

v co-ordinate of point 

Global v co-ordinate of an internal point 

v co-ordinate of Gauss point 

v co-ordinates of the boundary segment end points 

v co-ordinate relative to the origin for the particular integral 

y co-ordinate of the origin for the particular integral function 

y co-ordinate of the centre of curvature of a curved boundary segment 

y co-ordinates of the element end points 

y co-ordinate of first end point of current boundary segment 

y co-ordinate of an internal point 

Values of y direction integrals stored for domain integration 
y co-ordinate of last end point of current boundary segment 
y co-ordinate of the mid point between the ends of a curved segment 
y co-ordinates of the nodes 
y co-ordinate of point P 
Global co-ordinate of an internal point 
y co-ordinate of Gauss point Q 
y co-ordinates of the boundary segment end points 
y co-ordinate relative to the origin for the particular integral 
Intrinsic co-ordinate £; 

Values of the intrinsic co-ordinate $, at the Gauss points (normal quadrature) 
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5 Boundary Element Analysis of 
Plane Elastic Problems 


In this chapter a form of boundary element analysis for two-dimensional elastic stress analysis problems 
such as those outlined in Chapter 1 is presented. Only quadratic elements are considered. Three- 
dimensional problems are discussed briefly in Section 7.3. 

In Sections 1.2.5 and 1.2.6 it was shown that both plane strain and plane stress problems can be described 
in terms of a fourth-order biharmonic differential equation for Airys stress function. The relevant 
Equations are 1.71 and 1.77, both of which are special cases of Equation 1.85. Clearly, such problems are 
fundamentally different from potential problems, which are governed by the second-order Equation 1.84. 

Using an Airy stress function approach to solving plane elastic problems can be very appropriate 
when the boundary conditions are defined in terms of stresses. In more general problems of practical 
engineering interest, however, boundary conditions are typically defined in terms of a mixture of stresses 
and displacements, and a different approach is to be preferred. 
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Following the methodology developed in Chapter 2 for potential problems, what are required are a 
fundamental solution, equivalent to Equation 2.13, and a relationship equivalent to Greens symmetric 
identity, expressed in Equation 2.25. 

5.1 Fundamental Solution 

In the case of potential problems, the fundamental solution (Section 2.1) was in practical terms that due 
to a source of potential concentrated at a point in a solution domain of infinite extent in all directions. 
A mathematical requirement of a fundamental solution is that its value is singular (goes to infinity) at 
a point - the point where the source is located. 

In elasticity problems the corresponding fundamental solution is that due to a force concentrated at 
a point in an infinite domain, which again has the property of singularity at the point concerned. 
Immediately there is a substantial difference between the two classes of problems. The fundamental 
solution for potential problems involves only a scalar quantity: the source strength, and the effects of the 
point source in terms of potential distribution will be the same in all directions moving away from the 
point. In contrast, the fundamental solution for elasticity problems involves a vector quantity: the point 
force has both magnitude and direction. The effect of the force in terms of stresses and displacements 
is certainly not the same in all directions moving away from the point. 

In three dimensions, the fundamental solution is truly that due to a concentrated point force. In two 
dimensions, however, the problem is either one of plane stress, when the solution domain is very thin 
in the third dimension normal to the domain (Section 1.2.6), or plane strain, when the solution domain 
is very thick in the third dimension (Section 1.2.5). In either case, the fundamental solution needs to be 
thought of as that due to a force uniformly distributed along a line through the domain thickness in the 
third dimension. The force is a line force per unit thickness, which is applicable to either plane stress or 
plane strain. In the plane of the solution domain it still appears as a force acting at a point. 

Rather than deriving the fundamental solution from first principles, the approach adopted here is to state 
the result, and then show that this satisfies all the relevant requirements. Figure 5.1 shows a plane with 
both Cartesian (x,y) co-ordinates and polar ( r,0 ) co-ordinates, both with the same origin, and the 
angular co-ordinate measured anti-clockwise from the r direction. Stress components cr (radial direct 
stress), o ee (hoop direct stress) and a r6 (shear stress) are shown acting on a small region of the domain. 
As in Equation 1.1, shear stresses are complementary, so that o rQ — o Qr , and o t q is shown acting on 
all four faces of the region. 
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Figure 5.1 Stress components in polar co-ordinates 


Due to a line force of unit magnitude (per unit length) acting in the x direction at the origin, the three 
stress components at distance r from the origin and angular location 6 are given by 


(3+v*)cos 6 
Anr 


(5.1) 


a ee - 


(1—v*) cos 6 
Anr 


(5.2) 


a rd ~ 


(1—v*) sin 6 
Anr 


(5.3) 


The displacement components, u r in the radial direction and u e in the hoop direction, are given by 


u r 


(l+v*)(3—v*) cos 0 ^ /d\ 
AnE* \rj 


(5.4) 


u e = [(1 + v*) 2 - (1 + v*)(3 - v*)ln (£)] 


(5.5) 


where the length d is arbitrary, and remains to be chosen. It appears because the problem (of a 
concentrated force in an infinite medium) has had no displacement boundary conditions defined. 


The parameters E* and v* in these equations are the effective Youngs modulus and Poissons ratio for 
the elastic material of the domain. Their definitions depend on whether the problem is one of plane 
stress or of plane strain. For plane stress 


E* = E and v* = v 


(5.6) 


while for plane strain 


E* = 


E 

(1-v 2 ) 


and 


v 


V 

(1-v) 


(5.7) 
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In other words, under plane stress the effective properties are the actual properties, while under plane 
strain they are modified according to Equations 5.7. A plane strain problem can be treated as plane 
stress, provided the modified properties are used. 



Figure 5.2 Stresses acting on a circle about the point of force application 
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Figure 5.2 shows a circular region of radius r of the domain centred at the origin. The stresses acting 
on a small arc AB of the outer surface of the region of angular extent d 9 are o rr and o rQ as shown. 
For unit thickness of domain in the direction normal to the plane shown, the forces on AB in the radial 
and tangential directions are <j rr X rdO and a rd X rd(9, respectively, and the total force on the region 
in the negative x direction, which should be equal to the applied line force, is 

Force = / Q 27r (-( i rr cos 9 + a rQ sin0)rd0 (5.8) 

= J 0 2?r ^ [(3 + v*) cos 2 9 + (1 - v*) sin 2 9]d9 (5.9) 


Now 



(5.10) 


so that Equation 5.9 becomes 

Force = ^[(3 + v*)iz + (1 -v*)n] = ^(3 + 1) = 1 ( 5 . 11 ) 

The force is of unit magnitude, as required. An integral expression similar to Equation 5.8 can be used to 
show that the total force in the y direction is zero. Note that the requirement that the total forces on the 
circular region be constant, independent of radius, means that all the stresses are proportional to r -1 . 


In plane polar co-ordinates, the strains defined in terms of displacements (equivalent to Equations 1.2 
and 1.3) are 


e 


TT 


du r _ (l+v*)(3—v*) cos 6 

dr AnE*r 


(5.12) 


e ee 


U r 1 dug _ (1+v*) 2 COS 9 

r r dr AnE*r 


(5.13) 


_ 1 du r d _ 2(l+v*)(l—v*)sin 6 

^ r 0 r dQ r dr V r / AnE*r 


(5.14) 


Identical results should be obtained by applying the elastic constitutive equations (equivalent to 
Equations 1.20 to 1.25) to the stresses defined in Equations 5.1 to 5.3. 


Plane stress 

Under plane stress conditions and in the absence of temperature changes, the direct stress o zz normal 
to the plane of the domain is zero, and 

e r r =|(o rr ~ Vff g0 ) = ^ [-(3 + v) - v(l - v)] = - ^ ~ (5-15) 
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which, in view of Equations 5.6, is identical to Equation 5.12. Similarly 


e e e = \ O 00 ~ v <?rr) = ^ [(1 - v) + v(3 + v)] 


(1+v) 2 cos 9 
4nEr 


which is identical to Equation 5.13. Finally 



2(l+v) _ 2(l+v)(l—v)sin 9 

E 4nEr 


which is identical to Equation 5.14. 


(5.16) 


(5.17) 


Plane strain 

Under plane strain conditions, and again in the absence of temperature changes, the direct strain e zz 
normal to the plane of the domain is zero, from which 

e zz = \ Wzz ~ vOrr + oee )] = 0 (5.18) 

a zz = v(a rr + a ee ) (5.19) 

Then, using Equations 5.7 

1 1 

e rr = £ Wrr - v((Jee + O] = g [o'rr “ v(<J e 8 + VCTrr + V^gg)] 


= £ [(1 - V 2 )a rr - v(l + v)a 9e ] 


=?[ 


E* l ffrr 1 —v 099 


] =y t [°rr ~V*<Tgg] = 


(l+v*)(3—v*) cos 6 
4nE*r 


as in Equation 5.12. Also 


1 1 

&66 \.&99 ^(*7zz ^rr)] 1^99 v(v(J rr + VCTqq ^r)] 


= ^[(1 - v 2 )a de - v(l + v)a rr \ 


g, [&00 j _ y ^rr ] g* [&Q9 v °rr] 


(1+v*) 2 cos 9 
-1 

47r£’*r 


as in Equation 5.13. Finally 


_ a r g _ 2(l+v) _ 2(l+v*)(l-v*)sin 9 

e r9 ~ — ~ &r9 ~ 4nE*r 


as in Equation 5.14. 


(5.20) 


(5.21) 


(5.22) 
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Equations 5.15 to 5.17 and 5.20 to 5.22 show that the expressions for strains, under either plane stress or 
plane strain conditions, obtained from the stresses via the constitutive equations, are identical to those 
obtained from the displacements. With the stresses also satisfying equilibrium with the applied line force, 
this means that Equations 5.1 to 5.5 represent the true solution for a line force at a point in an infinite 
plane, in other words the fundamental solution for plane elastic problems. 

5.1.1 Displacement kernel functions 



These results now need to be generalised for line forces applied in both the x and y directions, to give 
the displacements and tractions in the same pair of directions. The term traction will be explained 
shortly. Figure 5.3 shows the displacements at a field point q produced by a line force of unit magnitude 
in the x direction at force point p. Displacements u r and u e are given by Equations 5.4 and 5.5. The 
corresponding displacement in the x direction at point q is 


u — u r cos 9 —u e sin 9 

2 p /W\ c!ki2 


(1 + v*)(3 - v*) cos z 9 (d\ si n z 9 

-In 


47r£* 


■©- 


4nE* L 


(1 + v*)(3 — v*) (d 
— In 


4nE* 


© 


(1 + v*) 2 - (1 + v*)(3 - v*)ln (-) 
(1 + v*) 2 sin 2 6 


4nE* 


(i+v*) 2 r(3-v*) 


4nE* L(l+v*) 


[§TS ta (7)- 1 + cos29 ] 


(5.23) 


It is convenient to choose the arbitrary length d such that 


u = 


(l+v*) 2 r (3-v*) 
4nE* l_(l+v*) 



+ cos 2 


(5.24) 


Now cos 9 = f x 


sin 9 = r y 


(5.25) 
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where f x and r y are the components in the % and y directions of the radius vector of unit length in the 
direction from point p to point q. Equation 5.24 can be written as 


U = U xx (p, q ) = 


(1+v*) 2 

[(3—v*) 

AnE* 

L(l+v*) 


In 


© 


+ V X T X 


r(p, q) =£ 0 


(5.26) 


U xx (p, q) is the displacement kernel function defining the displacement in the A' direction at q due to 
a unit line force in the x direction at p. 


Again from Figure 5.3, the displacement in the y direction at point p is 


v = u r sin 9 + u e sin 6 


(1 + v*)(3 — v*) cos 0 sin# /d\ cos 9 sin 9 


47 tE* 


© 


In - + 


4nE* 


v = U xy (p, q) = 


(i+v *) 2 

4ttE* 


cos 6 sin0 = 


( 1 +v *) 2 , . 
4 ? te* TxT y 


(1 + v*) 2 - (1 + V*)(3 - v*)ln 

(5.27) 


U X y (p, q ) is the displacement kernel function defining the displacement in the y direction at q due to 
a unit line force in the v direction at p. 
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For a line force of unit magnitude in the y direction at point p (Figure 5.3), similar results can be 
obtained. The angular position of field point q relative to the line of action of the force is now not 9 but 
— — 0^, so sin 0 becomes — cos 0, and cos 9 becomes sin 9. Using Equation 5.24, the displacement 

in the direction of the applied force at q is 


V = Uyy (p, q ) 


(i+v*) 2 r (3—v*) , 
4jtE* L(1+v») 



+ sin 2 


(l+v*) 2 p—v*) 
4 j iE* L(l+v») 



+ fyfy 


(5.28) 


r(p, q) * 0 


Using Equation 5.27, the displacement in the direction at 7r/2 anti-clockwise from the direction of the 
applied force at q is 

( 1+V *) 2 • n 

—u =-sin 9 cos 9 

4nE* 

and u = U yx (p, q) = r y f x (5.29) 

The four equations, Equations 5.26 to 5.29, define the four displacement kernels. They can be expressed 
elegantly in a single equation if tensor notation is introduced. For present purposes, however, it is more 
convenient to have them in their more explicit, if more lengthy, forms. The first subscript of U indicates 
the direction of the unit line force at p, while the second subscript indicates the direction of the resulting 
displacement at q. 

5.1.2 Traction kernel functions 

The displacement kernel functions are expressed in terms of displacements at the field point in the global 
co-ordinate directions. Similar results are required for stresses. For this purpose, the concept of tractions is 
introduced. A traction is the resultant stress, or force per unit area, on a surface in a particular direction. 



Figure 5.4 Tractions and stresses at a surface within a domain 
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Figure 5.4 shows both the Cartesian stress components and the equivalent tractions, t x and t y in the * and 
y co-ordinate directions, acting on a surface of a small triangular piece of material within a domain. The 
unit outward normal vector to the surface, n, is inclined at angle 8 to the v direction. For equilibrium 
of forces in the x direction (bearing in mind that the domain is of unit thickness) 

t x x (AC) = a xx x (AB) + a xy x (BC) 

t x = a xx cos S + a xy sin 5 = a xx n x + a xy n y (5.30) 

where n x and n y are the components in the co-ordinate directions of the unit normal vector, and AC, 
AB and BC are the lengths of the sides of the small triangle. Similarly, for equilibrium of forces in the 
y direction 

ty X (AC) = (Jyy X (BC) + X (AB) 

ty = (Jyy Sm 8 O X y COS 8 = (Jyy Ti y (J X yYl X (5.3 1 ) 



Figure 5.5 Stress components in polar co-ordinates 

Figure 5.5 shows the polar co-ordinate stress components at field point q due to a line force in the x 
direction at p. Also the equivalent tractions in the global co-ordinate directions,^ and t y , acting on a 
surface of a small triangular piece of domain whose outward normal is inclined at angle y to radius r, and 
angle 8 to the v direction. For equilibrium of forces on the right-angled triangular piece in the v direction 

t x x (AC) = <j rr x (AB) cos 6 - o rQ x (AB) sin 9 - a ee x (BC) sin 9 + a r6 x (BC) cos 9 

t x = o rr cos y cos 9 — o rQ cos y sin 9 — a ee sin y sin 9 + a rQ sin y cos 9 
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With the stresses defined by Equations 5.1 to 5.3, this becomes 

(3 + v*) cos 2 9 cos y (1 - v*) sin 2 9 cos y 
x 4nr 4nr 

1 

= - -— [(1 - v*) + 2(1 + v*) cos 2 0 ] cos y 
4nr 

tx = - ^ [(1 - v*) + 2(1 + v*) cos 2 0] £ (5.32) 

where — = cos y. The traction kernel function for defining the traction in the % direction at q due to 

dn 1 

a unit line force in the x direction at p is 

1 c\y 

T xx (p,q) = ~ —[(l~v*) +2(1+v*)r x r x ] — r(p,q) + 0 (5.33) 

Similarly, for equilibrium of forces in the y direction 


t y x (AC) = o rr x (AB) sin 6 + o r6 x (AB) cos 0 + a ee x (BC) cos 0 + a rd x (BC) sin 6 
t y = a rr cos y sin 6 + o r6 cos y cos 6 + a ee sin y cos 9 + o r6 sin y sin 9 

= —— [(3 + v*) — (1 — v*)] cos 9 sin 9 cosy + ^ v - siny (5 34 ) 

4nr v ' 1 4 nr ‘ v ' 
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Now sin y — sin(5 — 0) = sin S cos 9 — cos S sin 9 

— n y f x — n x f y (5.35) 

where n x and n y are the components in the % and y directions of the unit outward normal to surface 
AC at point q. The kernel function for defining the traction in the y direction at q due to a unit line 
force in the % direction at p is 

T X y (p, <?) = - [:2(1 + v*)f x f y £ - (1 - v*)(f x n y - r y n x )] (5.36) 

r(p, q) =£ 0 

For a line force of unit magnitude in the y direction at point p (Figure 5.5), similar results can be 
obtained. The angular position of field point q relative to the line of action of the force is now not 9 but 
— — 0 ), so sin 9 becomes — cos 0> and cos 9 becomes sin 9. Also, n x becomes n y and n y becomes 

—n x . Using Equation 5.32, the traction in the direction of the applied force at q is 

h = ~ ^ K 1 - v *) + 2(1 + v*) sin 2 #] £ (5.37) 

and T yy (p,q) = [(1 - v*) + 2(1 + v*)f y f y ]^ r(p,q) * 0 (5.38) 

Using Equation 5.34, the traction in the direction at n/2 anti-clockwise from the direction of the 
applied force at q is 

1 r dr 

—t x = — -— 2(1 + v*) (—sin0) cos0 -(1 — v*)(— n x sin# — n v (— cos#)) 

47rr 1 dn 7 

t x = — [2(1 + v*) sin # cos # ^ — (1 — v*)(n x sin # — n y cos 0)J (5.39) 

and ^(p,?) =-^[2(1+ v*)f y f x ^-(1-v*)(r y n x -f x n y )] (5.40) 

r(p, q) * 0 

Again Equations 5.33, 5.36, 5.38 and 5.40 for the traction kernel functions can be generalised in a single 
formula using tensor notation. 

The displacement and traction kernel functions defined here for elastic stress analysis are equivalent to 
those defined in Equations 2.13 and 2.34 for potential problems. There is a considerable difference in 
complexity, reflecting the vector nature of a concentrated force against the scalar nature of a point source. 
But, there are important similarities in that both pairs are functions of either In (r -1 ) (displacements 
and potential) or r _1 (tractions and potential gradient). The techniques of dealing with such functions 
developed for potential problems will be applicable to elastic problems. 
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5.2 Boundary Integral Equations 

The relationship equivalent to Greens symmetric identity for potential problems required to develop the 
boundary integral equation for elastic stress analysis problems is Betti s reciprocal theorem. Suppose that 
an elastic body is subject to a number of forces, F k , at particular points and in particular directions on 
its surface, and that the corresponding displacements at the same points and in the same directions as 
the forces are u k . Suppose also that the same body is separately subject to another independent set of 
forces, F/*, at the same points and in the same directions, producing corresponding displacements u\. 
Bettis theorem states that 

F k u k = Zfc FkU k (5.41) 

It is a form of virtual work principle: the total virtual work done by the first set of forces moving through 
the second set of displacements is equal to the virtual work done by the second set of forces moving 
through the first set of displacements. 

For present purposes, the first set of forces and displacements can be those associated with the particular 
problem to be solved, and the second set with the fundamental solution. Also, a force F, which must in 
reality act over a finite area, is equivalent to a traction, t, and the summation in Equation 5.41 becomes 
an integral over the boundary surface, S 


f s t k u* k dS = { s t k u k dS (5.42) 

Suppose that some point within the solution domain serves as the origin and point of application of 
the concentrated forces (in the two co-ordinate directions) for the fundamental solution, as shown in 
Figure 5.6 



Figure 5.6 A solution domain, including a small region of radius surrounding the point p 

Now, the fundamental solution is not valid at the point p itself. Consider therefore a small internal 
boundary, S f , of radius s centred at p. Equation 5.42 can be written for the region between S and S £ , 
which excludes p where both displacements and tractions are singular 

f s t k u* k dS + f Se t k u* k dS = f s t* k u k dS + f Se t* k u k dS (5.43) 
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The next step is to consider what happens as the radius £ shrinks to zero. On the circle, a displacement 
kernel u k (Equations 5.26 to 5.29) takes the general form of 


u k 


k=A(0)ln©+/ 2 (0) 


(5.44) 


where /i(0) is zero for U xy and U yx . If £ is small the value of traction t k is effectively constant over the 
boundary and equal to its value at point p. Therefore 


f t k u* k dS = t k (p) f /i(0)ln(-)£d0 + t k (p) f f 2 (d)add 
Js £ Jo V£/ J 0 

= tfc(p) £ In (I) Ci + t k (p)e C 2 


where C[ and C 2 are constants. The second term obviously goes to zero as a -* 0. Similarly, 
lim £ ^ 0 [gin 0)] = 0, and consequently 


f t k u* k dS = 0 


Also, a traction kernel tt must be such that 


(5.45) 


Jc t* k dS = l 


(5.46) 
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if the traction is in the direction of the unit concentrated force at p (and zero if it is at right angles to 
the force). Therefore, if the force is in the % direction 

f s t k udS = u(p) (5.47) 

d £ 

and as the boundary S £ shrinks to zero size at the point p, Equation 5.43 becomes 

u(p) = f s t k u* k ds - f s t* k u k cLS (5.48) 

The concentrated force in the v direction at p creates displacements and tractions in both the v and y 
directions at all points on the boundary of the domain, as defined by the kernel functions. Therefore 

<v) = [ [t x (QW xx (p,Q) + t y (Q)U X y(p,Q)]dS(Q) 

J S 

- f s [u(Q)T xx ( p, Q ) + v(Q)T xy (p, <3)]dS(Q) (5.49) 

Similarly, if a unit line force is applied in the y direction at p 

v(p ) = [ [t x (Q)U yx (p, Q) + t y (Q)U yy (p, Q)]dS((?) 

■>s 

~ f s [u(Q)T yx (p, Q ) + i KQ)T yy (p, Q)]dS(Q ) (5.50) 

The point Q is a field point on the boundary which moves along the boundary as the integrations proceed. 

The displacement and traction kernel functions are known, so provided the values of both the 
displacements and tractions are known at every point along the boundary, Equations 5.49 and 5.50 
provide a means of calculating the solution displacements at any point within the solution domain. 
They can be differentiated to find strains at p, and hence stresses via the elastic constitutive equations. 

In engineering practice, stresses and strains at points within an elastic solution domain are rarely required. 
In all but a few classes of problems the largest individual stress components and the most intense states 
of stress are located on a boundary. So for present purposes, data at internal points are not computed. 

Let point p be taken to the boundary at a point P, and again consider a small region of exclusion with 
boundary S £ of radius centred at P. While Equation 5.45 still holds, Equation 5.46 becomes 

f t%. dS = constant (5.51) 

*->£ 
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if the traction is in the direction of the concentrated force at P, where the constant is no longer unity. 
Indeed, if the boundary at P is smooth, the value of the constant is V2, because exactly half of the 
complete circle is within the solution domain. If the boundary is not smooth, however, the value of the 
constant is difficult to evaluate directly. Similarly, if the traction is in the direction at right angles to the 
direction of the concentrated force at P, Equation 5.51 also applies, where the constant is no longer 
zero. Equations 5.49 and 5.50 become 

C xx (P)u(P) + C xy (P)v(P) = J [t x ( Q)U XX (P, (?) + ty ( Q)U xy (P, <?)]dS(<?) 

- I s [u(Q)T xx (P, (?) + v(Q)T xy (P, <?)]dS(<?) (5.52) 

C yx (P)u(P) + C yy (P)v(P) = [ [t x (Q)U yx (P,Q) + t y (Q)U yy (P, <?)]dS((?) 

Js 

- Is [u(Q)T yx (P, (?) + v{Q)T yy (P, <?)]dS(<?) (5.53) 

where C xx (P),C xy (P), C yx (P) and C yy (P) are constants. The free terms containing these constants 
contribute only to the coefficients at the diagonal of the relevant matrix in the numerical implementation 
of the method, and these coefficients can always be evaluated indirectly. 

Equations 5.52 and 5.53 are now a pair of boundary integral equations for point P as the point at which 
line forces of the fundamental solution are applied, the first for force in the x direction and the second 
for force in the y direction. They can be used to determine the distributions of the displacements and 
tractions over the boundary. 

5.3 Discretisation of the Boundary Integral Equations 

In general, Equations 5.52 and 5.53 cannot be solved analytically, and some form of numerical method 
must be employed. Boundary displacements and tractions are found not as a continuous algebraic 
functions of position along the boundary, but as numerical values at a finite number of discrete points 
on the boundary. The boundary may be subdivided into small pieces or boundary elements. Associated 
with each element are one or more of these points: the nodes or nodal points. The distributions of 
displacements and tractions over the elements are defined in terms of nodal point values by suitable 
interpolation functions. 
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Boundary integral Equations 5.52 and 5.53 are applied to each of the N nodal points P in turn, in the 
discretised form 


M 

c xx (P)u(P) + C xy (P)v(P) + £J [iu(Q)T xx (P, (?) + v(Q)T xy (P, <?)]dS«?) 


?71 = 1 


M 

-XL [tx(Q)U XX (P, (?) + ty ( Q)U xy (P, (?)] dS(<?) 


m= 1 JSm 


(5.54) 


M 

C yx (P)u(P) + C yy (P)v(P) + ^J [u(Q)T yx ( P, (?) + v(Q)T yy (P, <?)]dS(<?) 


771=1 


M 

-XL [t x (Q)U yx (P,Q) + t y ((?)t/ yy (P,<?)] dS((?) 


m=l ' /5m 


(5.55) 
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The total number of boundary elements is M, m is an element counter, and S m is the piece of the boundary 
occupied by element number m. The details of how integration is carried out over an individual element 
depend on the type of element involved. Equations 5.54 and 5.55 represent a set of 2N linear equations, 
where N is the number of nodal points (N = M for constant or linear boundary elements, N = 2M for 
quadratic elements). It is convenient to arrange them in the order of the node numbers, with for each 
node the equation corresponding to the v direction force first, followed by that for the y direction force 

[A][u] = [B][t] (5.56) 

Matrices [A] and [B] are square and contain known constant coefficients, in general all of which will 
be non-zero. The column vectors [ u] and [t ] contain the nodal point values of displacements and 
tractions, two components for each at each node. At each node in each direction either the displacement 
or traction will be unknown and the other known, or there will be a linear relationship between them. 
The equations can be rearranged, taking all unknown quantities to the left hand side, known quantities 
to the right hand side, giving a set of linear equations in a familiar form 

[A*][x] = [B*][y] = [b] (5.57) 

where [A*] and [/>* ] are modified coefficient matrices and [fi] is a column vector of known coefficients. 
This set can be solved for the 2N unknowns v at the nodes, meaning that the displacements and tractions 
are known at every nodal point on the boundary. Given this information, stresses at the nodes can also 
be found. 

5.4 Boundary Conditions and Surface Stresses 

The most straightforward type of boundary condition met in practice is where tractions are prescribed. 
In other words, for a particular node the traction components in both of the global co-ordinate directions 
are known. These values contribute to the right hand side of Equations 5.57, the solution of which gives 
the displacement components at the node. 

Traction boundary conditions are most likely to be prescribed in the form of stresses on the boundary, 
in particular the direct stress normal to the boundary and the shear stress along it. 
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Figure 5.7 Stresses and tractions at a point on the boundary 

Figure 5.7 shows the normal stress a nn and shear stress a sn acting at a point on the boundary. Note 
the sign convention for shear stresses: a positive shear stress acts in the positive Sdirection, along the 
boundary keeping the domain to the left. The outward unit normal vector n at the point concerned is 
inclined at angle y to the v global co-ordinate direction, and has components n x and fi y in the v and 
y directions. Because both stresses and both tractions (resultant stresses) act on the same surface, they 
can be resolved like forces 

t x Onn COSy G sn siny ^nn^x &sn^y (5.58) 

t y = a nn siny + a sn cosy = a nn u y + o sn n x (5.59) 

With displacement boundary conditions, for a particular node the displacement components in both 
of the global co-ordinate directions are known. These values must be taken from the left hand to the 
right hand side of Equations 5.57 to contribute to [b\. The solution of the equations gives the traction 
components at the node. 

It should be noted that tractions are often discontinuous at a corner, and indeed at any point where either 
the direction of the outward normal to the boundary or the magnitudes of the stresses, change abruptly. 
Figure 5.8 shows a right-angled corner of a domain with an applied normal stress on one side of the 
corner, and stress free on the other side, a situation that occurs commonly in practice. Label A indicates 
the point of the corner. Consider the two points B 1 and B 2 on the two sides of A. At B x the tractions are 
t x = 0 ,t y = 0, whereas at B 2 they are t x = a 0 , t y = 0. These traction values remain unchanged as B l 
and B 2 approach A. Clearly, there is a discontinuity in traction t x at A, which raises the question as to 
what is meant by the tractions at such a point. 
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Gnn ~ &sn — 0 A 



Figure 5.8 Tractions at a corner 

In practice, there is no difficulty if stresses are prescribed on both sides of a corner, provided the 
integrals which must be evaluated over both sides of the corner employ the relevant traction values, 
without reference to any (undefined) corner value. The displacement components at the corner, which 
are not discontinuous, are the unknowns to be found. The case of displacements prescribed on both 
sides of a corner is not of great practical interest, corresponding to a form of rigid body displacement. 
At least for uniform displacements the tractions would be zero. If one side of a corner has displacements 
prescribed, and the other stresses, then it is convenient to regard the corner as being part of the prescribed 
displacement side, and treat the unknown tractions as the tractions applicable to this side (the tractions 
on the other side being known). 
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Other types of boundary conditions create greater complexity For example, a line of symmetry has fixed 
(usually zero) displacement normal to the line, zero traction along it. Such a mixture of displacements 
and tractions can be dealt with, but corners where a line of symmetry meets a part of the boundary with 
either prescribed displacements, prescribed stresses or indeed another line of symmetry, require special 
care. Further possible types of mixed boundary condition include, for example, flexible supports where 
tractions are related to displacements. 

It is important to note that in boundary element methods for elastic stress analysis problems it is not 
possible to apply point force boundary conditions, which would give rise to infinite stresses. This is in 
contrast to other numerical methods such as finite elements. It does, however, reflect physical reality 
in which even concentrated forces on a body are applied over small but finite areas. Similarly, point 
displacement constraints are not permitted if they would require point forces to maintain them. On the 
other hand, such constraints applied to a domain which is already in force equilibrium, simply in order 
to prevent movement of the domain as a rigid body, are permitted. Indeed, they are necessary to prevent 
such a problem giving rise to a singular [A*] matrix in Equations 5.57. This can be a very useful facility, 
and will be demonstrated in the problems considered in Chapter 6. 

Point constraints are straightforward to apply. If, say, a particular node needs to be constrained in the v 
direction, the first of the two equations corresponding to that node is modified. All the coefficients in 
that equation, including that in the right hand side vector, but excluding the coefficient on the diagonal 
of matrix [A*], are set to zero. The diagonal coefficient is set to one. 

The primary results of the analysis are the displacements and tractions at the nodes. Stresses are likely 
to be of much greater interest in practice than tractions. Referring to Figure 5.7, the normal and shear 
stresses can be obtained from the tractions as 


o nn — t x cosy + ty siny — t x 7i x T ty^iy 


(5.60) 


°sn = ~t x siny + t y cosy = -t x n y + t y n x 


(5.61) 


The other stress component of interest is the direct stress in the direction along the surface, g S s- This 
cannot be found from the tractions, but can be found from the displacements. If u s is the displacement 
along the boundary in the S direction, the direct strain there is given by 


p — —- 

ss dS 


du s 


(5.62) 


But this strain can also be found from the stresses and elastic properties 


e ss £.* (&SS V ) 


from which 



(5.63) 
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Given the three stress components, cT nn ,a ss and a sn , other stress parameters can also be found, such as 
the principal stresses at the point. Often of greater practical use is a single stress, reflecting the overall 
intensity of the state of stress and hence the likelihood of yielding or failure. The choice of this single or 
equivalent stress therefore depends on the criterion adopted for yielding or failure. At least for ductile 
materials, the most commonly used is the von Mises equivalent stress, which in the present situation 
is defined as 

C e V ^nri "h 0"ss Onn Tss "h 3 0" sn (5.64) 

To test for yield or failure, this equivalent stress is compared with the yield or failure stress for the 
material concerned. 

5.5 Quadratic Boundary Elements 

Experience of potential problems in Chapters 3 and 4 showed that, although constant boundary elements 
are easier to implement, quadratic elements give much more accurate results for the same fineness of 
discretisation. The only exceptions to this were when results were required at points within a solution 
domain close to a boundary, or when the domain was very narrow so that opposite sides of the boundary 
were close together. For elastic stress analysis problems attention is confined here to quadratic boundary 
elements. This is particularly because such problems are often concerned with rapid stress variations in 
local regions such as stress concentrations or cracks. Results at internal points are rarely required, but 
care will need to be taken with long slender domains. 


P 



Figure 5.9 Quadratic element discretisation of a two-dimensional solution domain boundary 
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Figure 5.9 shows a typical arrangement of quadratic line elements on parts of the boundary of a two- 
dimensional solution domain. Each element has three nodes, one at each end and one at its centre. While 
the end nodes are shown as solid circles, those at the centres of elements are shown as open circles. The 
relevant boundary integral equations are Equations 5.54 and 5.55. The total number of nodal points 
is twice the total number of elements, and there are two equations per node, hence 4 M equations are 
generated, M being the number of elements. The integrations involved in the evaluation of the equation 
coefficients must in general be carried out numerically. 
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Figure 5.10 shows a typical curved element, exactly as used for potential problems. Its three nodes are 
numbered in the direction of integration from 1 to 3, these numbers being local to the particular element. 
The second node is at the centre of the element as measured along its curved length. A local intrinsic 
co-ordinate, <f, follows the curved shape of the element. It has its origin at the centre node: that is, = 0 
at node 2, and takes the values -1 and +1 at nodes 1 and 3, respectively. 

The distributions of displacements and tractions are approximated as quadratic functions of position 
along the element in terms of values at the three nodes. For example 


u(0 = Wi(0«i + N 2 (0«2 + N 3 (0«3 = 2=i^c(0“c (5.65) 

Shape functions are as defined in Equations 2.50, 2.51 and 2.52 

= (5.66) 

N 2 (0 = 1-^ (5-67) 

N 3 (.0 = \ttf + 1) (5.68) 

Variations of the global co-ordinates are similarly defined 

xCO = Ni (0*1 + N 2 ( 0*2 + N 3 (0*3 = Z C 3 =1 N c (0*c (5.69) 

y (0 = n , (Oyi + n 2 (Oy 2 + n 3 (0 y 3 = S c 3 =i n c (O y c (5.70) 


The rates of change of the global co-ordinates with respect to <f are 


dx _ y3 d N c tf) dy_y 3 d NdO 

df Zjc=1 df c? df Z,c=1 df 


and the derivatives of the shape functions are 


dVi(£) _1 d V 2 (O diVsCO r , 1 

df ^ 2 ' df ^ 7 df ^2 


As in Equation 2.59 the vector normal to the boundary can be found as 


(5.71) 


(5.72) 


n = n x i + 

where n x and n y are 


c 1 a i dy . dx . 

"W=SAlc = -.-- 7 

the components in the v and y directions of the vector normal n. 


(5.73) 
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In order to integrate round the boundary it is necessary to change from global co-ordinate to local 
intrinsic co-ordinate £ within each element, by means of the Jacobian of transformation 

-1 - + <5 - 74) 

which can also be used to define the components of the unit normal 


n = 






Tlx * I n y * » . yv . 

= m l+ m' = n * l + n y'' 


(5.75) 


Introducing the transformation from global to intrinsic co-ordinate, the boundary integral equations, 
Equations 5.54 and 5.55, become 


c xx (P)u(P) + C xy (P)v(P) + Z / [u(Q)T xx (P, (?) + v(Q)T xy (P, (?)] /(Odf 

m =1 — ^ 

M +1 

-XL [t x (Q)U xx (P,Q) + t y (Q)U xy (P,Q)]j(Odf 

m= 1 — ^ 

(5.76) 


0* (P)u(P) + C yy (P)v(P) + ]T J [u(Q)T yx ( P , (?) + v(Q)T yy (P, (?)] /(Odf 

m=l — 1 

M x 

-XL [* t x (Q)U yx (P,Q ) + t y ((?)i/ yy (P, (?)] 7(Odf ( 

m=l _1 

Then introducing the parametric representation of Equation 5.65 for the displacement and traction 
distributions 


M 3 


Cxx (P)u(P) + C xy (P)V(P) + X Z / (P ’ Q) + ^ (P» <2)R (0 /(Odf 

m=l c=l _1 

M 3 +1 

-ZZ L [it x }cU xx (p, (?) + {t y } c 0y (p, (?)] tf c (0/(0d? 

m=lc=l _1 


(5.78) 


M 3 


C yX (P)li(P) + Cyy (P)V(P) + Z j / ^ + ^ ^ 7(0df 

771 — 1 C — 1 _1 

M 3 +1 

= ZZ (P. 0 + {ty}Uyy (P. «] 


771—1 C = 1 


(5.79) 
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where c is the number (1, 2 or 3) of the node in element number m. The first (traction) kernel functions 
on the left hand sides of the equations and the second (displacement) kernels on the right hand sides 
are given by Equations 5.33, 5.36, 5.38, 5.40 and 5.26 to 5.29. 

5.5.1 Points P and Q not in the same element 

The required integrations are performed using Gaussian quadrature, described in Appendix A. Provided 
point P is not in the element over which integration is required the process is straightforward. 

5.5.2 Points P and Q in the same element 

The case where P is in the relevant element requires some care, because when and are coincident the 
kernel functions are singular. The orders of the singularities are r -1 and ln(r -1 ). Provided that P is not 
the cth node of the element, however, the shape function N c (<f) goes to zero at P, and the products of 
kernels and shape functions are not singular at P, and may be integrated normally. 

If P is the cth node of the element then the terms at the diagonal of matrix [A], involving both the first 
(traction) kernel functions and the free term constants C xx (P),C xy (P),C yx (P), and C yy (P), are difficult 
to evaluate directly. They can, however, be evaluated indirectly. If there are no tractions applied to the 
boundary of a domain, Equations 5.56 become 


[A][u] = [0] 


(5.80) 
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Possible outcomes of this lack of loading include rigid body displacement of the entire domain in the x 
direction by a unit distance (all x-direction diplacements take the value 1, all y-direction displacements 
take the value 0). Also, rigid body displacement of the entire domain in the y direction by a unit 
distance (all y-direction displacements take the value 1, all x-direction displacements take the value 0). 
This means that the sum of all the coefficients multiplying x-direction displacements in any row of [A] 
should be zero, and, separately, the sum of all the coefficients multiplying y-direction displacements in 
any row of [A] should also be zero. These conditions allow the coefficients of displacements at point P 
(four coefficients in total, two in each of the two equations) to be determined by summations equivalent 
to Equations 2.47 and 2.66. 

Still considering P to be the cth node of the element in Equations 5.78 and 5.79, the second (displacement) 
kernels require special treatment to deal with the ln(r -1 ) singularity, using a modified form of Gaussian 
quadrature. The procedures follow closely those described in Section 2.5.2, but are repeated here for 
convenience. The radial distance from P to a point Q at (x, y) can be arranged in the form 

r(P, Q) = r]R(t;) (5.81) 

where i?(0 is a known function and 77 is a modified intrinsic co-ordinate with its origin at P. The form 
of R (O depends on which of the three element nodes P is located at. Co-ordinate 77 is chosen in each case 
to conform to the requirements of Gaussian quadrature involving a logarithmic function (Appendix A). 

P at the first node 

If P is at the first node of the element 

r{P,Q) = [(x-xi) 2 + (y-yi) 2 r ( 5 - 82 ) 

and x and y are defined in terms of intrinsic co-ordinate <f by Equations 5.69 and 5.70. So 

X - Xi = [N t (0 - l]Xi + N 2 (0*2 + N 3 (0*3 (5.83) 

y-yi = [tfi(0 - i]yi + n 2 (Ov 2 + w 3 (Oy3 (5.84) 
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Intrinsic co-ordinate 77 is chosen to range from 0 at the first node (<f = — 1) to 1 at the third node 
(<f = +1), so that 


V = ^ + D ^=2 

(5.85) 

[NdO - 1] = | (f 2 - ? - 2) = | (f + l)tf - 2) = J/tf - 2) 

(5.86) 

iv 2 (0 = 1 — f 2 = 277(1 — f) 

(5.87) 

= + 0 = »rf 

(5.88) 


Therefore 


X-X 1 = r][^ - 2)x 1 + 2(1 - 0*2 + f*3] (5-89) 

y-yi = vl(f - 2 )yi + 2(1 - 0y2 + M ( 5 - 9 °) 

«(0 = {[(? - 2)x x + 2(1 - 0*2 + f* 3 ] 2 + [«■ - 2)y x + 2(1 - Oy 2 + fy 3 ] 2 }* ( 5 - 91 ) 


The integrals of the second kernels in Equations 5.78 and 5.79 can be expressed in the form 
'+1 r+f 


f u xx (p,q)n c (OJ( 0^ = f 


(l + v*) 2 

[(3-0 

4nE* 

(1 + V*) 


-In 


(r(P, Q)} 


+ r x r x 


N c (OJ(Od( 


KiK 2 J In Q N c (OJ(0 |^| drj + K L K 2 f ^ In N c (OKO df + K, f f x f x N c (OJ(Od{ (5.92) 


f U xy (P,Q)N c (OJ(0 df = J + 


(1 + v ) 2 „ 

-r r r,, /V,. 

4nE* x y c 


r+l 

(07(0df =K t J ^ r x r y N c (OKOdZ (5.93) 


J + U yx (P,Q)N c (OJ(0d( = j + r y r x N c (OJ(0df =K t J + f y f x N c (07(0^(5.94) 


f Uyy(P,Q)N c (Omdf= f 


r +1 (1 + v*) 2 

P‘ v ’V 1 

J_! 47 tE* 

(1 + v*) " Vr(P, (2) 

r +1 , 

1 \ r 

+ K x K 2 J ^ In [— J N c (OJ(0 df + Ki J 


+ Tyfy 


NciOKO 


where O = ( ~ 1+v ^ and ^2 = ^ v ^ 

1 4jr£* Z (l+v*) 


(5.96) 


are constants. 


Because R (O is not zero within the range of integration, all of the integrals except those involving In 
(-J can be evaluated by normal Gaussian quadrature. Those involving the singular logarithmic function 
In y-j can be evaluated using the appropriate quadrature formula described in Appendix A. 
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P at the second node 

If P is at the second node of the element 

(5.97) 

(5.98) 

(5.99) 


r(P, Q ) = [(x - x 2 ) 2 + (y- y 2 ) 2 F 
x - x 2 = N.COxi + [N 2 (0 - l]x 2 + N 3 (0*3 
y-y 2 = N^Oyi + [n 2 (0 - i]y2 + N 3 (Oy3 


For integration purposes, the element needs to be divided into two regions, from the second node to 
the third node (<f = 0 to 1) and from the second node to the first node (<f = Oto —1). Between the 
second and third nodes the intrinsic co-ordinate is chosen as 





(5.100) 


antUViCO = — D = ~ 1) 


(5.101) 


[N 2 (0 -1] = -<r 2 = -mS 


(5.102) 


N 3 (0=\ttf + D=^i(f + D 


(5.103) 
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Therefore 


x - X 2 = 771 - 1 )*! - fx 2 + i(<f + 1 )* 3 ] 

y - y2 = Vi [j (? - i)yi - £y 2 + ^ (f + to] 


fl(0 = {[^ (f -to - ^ (f + to] + (f - to - ?y 2 + ^ (f + to 

Between the second and first nodes the intrinsic co-ordinate ;/ 2 is chosen as 

d? 


P2 = 


*?2 


= 1 


and P(0 is as defined in Equation 5.106. 


1 



(5.104) 

(5.105) 

(5.106) 


(5.107) 


The integrals of the second kernels U xy (P, Q) and U yx (P, Q ) in Equations 5.78 and 5.79 are as in 
Equations 5.93 and 5.94. The integrals of the other two can be expressed in the form 


0 


U xx (P,Q)N c (OKOdf 
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The third and fourth integrals on the right hand sides of Equations 5.108 and 5.109 can be evaluated by 
normal Gaussian quadrature, while the first and second involve the singular logarithmic function and 
must be evaluated using the appropriate quadrature formula. 
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P at the third node 

If P is at the third node of the element 

r(P,Q)= [(x - x 3 ) 2 + (y - y 3 ) 2 ]2 (5.110) 

x - x 3 = N^Oxi + N 2 (Ox 2 + [N 3 (0 - l]x 3 (5.111) 

y - y 3 = Ni(Oyi + N 2 (OV 2 + [N3OO - l]y3 (5.112) 

Intrinsic co-ordinate r\ is chosen to range from 0 at the third node (<f = +1) to 1 at the first node 
(<f = —1), so that 


n = \0-0 %=* 

(5.113) 


(5.114) 

N 2 (0 = l-f 2 = 2r?(l + 0 

(5.115) 

[N 3 (0 ~ 1] = ^(f 2 + f - 2) = i(f - l)(f + 2) = -j/tf + 2) 

(5.116) 


Therefore 


X - X 3 = T][-. fx x + 2(1 + Ox 2 - (f + 2 )x 3 ] (5.117) 

y - y3 = pHti + 2(1 + Oy2 - (? + 2)y 3 ] (5.118) 

R(0 = {[-fXi + 2(1 + Ox 2 - (f + 2 )x 3 ] 2 + [-f yi + 2(1 + Oy 2 - (f + 2)y 3 ] 2 F ( 5 - 119 ) 


The integrals of the second kernels U xy (P, Q) and U yx (P, Q ) in Equations 5.78 and 5.79 are as in 
Equations 5.93 and 5.94. The integrals of the other two can be expressed in the form 
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The second and third integrals on the right hand sides can be evaluated by normal Gaussian quadrature, 
while the first involves the singular logarithmic function and must be evaluated using the appropriate 
quadrature formula. 


5.6 Scaling 

As in the case of two-dimensional potential problems (Section 2.6), it is necessary to scale all distances 
between nodes to be less than unity. This is done by dividing by the maximum dimension of the problem 
(the maximum distance between any pair of nodes). Displacements have the units of length, so these 
are also scaled by dividing by the maximum dimension. Similarly, stresses and tractions are scaled by 
dividing by the value of Youngs modulus. Once the solution to the scaled problem has been obtained, 
the scaling is removed. 


5.7 Solving the Linear Equations 

As for potential problems (Section 2.7), the linear equations arising from the boundary element analysis 
are most appropriately solved by direct methods such as Gaussian elimination. The method is described 
in detail, including an appropriate computer subprogram, in Appendix B. 


5.8 Body Forces 

Many potential problems involve the Poisson form of governing differential equation rather than the 
simpler Laplace form. As a result, a particular integral has to be added, which can be done by modifying 
the boundary conditions for the Laplace problem, as explained in Section 2.8. In elastic stress analysis 
problems the equivalent situation arises if there are significant body forces. The most common of these 
is gravity, but other examples include centrifugal loading. Gravitational effects only become significant 
in physically very large components and structures, such as bridges and dams. In the large majority of 
problems body forces can be neglected, and they are not considered in detail here. 


Problems 

5.1 Under what conditions are the states of plane stress and plane strain indistinguishable? Is this 
of any practical significance? 

5.2 Find the particular forms of the typical displacement and traction kernel functions, 

and T xy (p,q ) for an incompressible material under plane 

strain conditions. 
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5.3 


The typical point on a boundary shown in Figure 5.7 is subject to flexible boundary conditions 


k n Un G sn k s U s 


where u n and u s are the displacements in the n and S directions, and k n and k s are normal 
and shear stiffnesses. Define the boundary conditions for the displacements and tractions in the 
global co-ordinate directions. 


5.4 


Assuming that gravity acts in the negative y direction, show that the following distributions of 
stresses and displacements provide a possible particular integral for this body force 

°yy ~ pay > 

u — 0 , 

where p is the material density, and g is the acceleration due to gravity. 


°xx — v °yy 


°xy = 0 


v — 


pgy 

2 E* 


(l-v*0 


5.5 If a body whose material has a density of p is rotated at an angular velocity of co about the y 
axis, the centrifugal body force per unit volume generated is pco 2 x in the x direction. Show that 
the following distributions of stresses and displacements provide a possible particular integral 
for this body force 
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6 Quadratic Boundary 

Element Program for Plane 
Elastic Problems 


In this chapter a Fortran computer program to implement the quadratic boundary element formulation 
for two-dimensional elastic stress analysis problems developed in Chapter 5 is presented and described 
in detail. It is then used to solve a range of problems to demonstrate the capabilities of the method. 

For readers who prefer to use Matlab, a translation is provided in Appendix E. 

As far as possible the program structure, file names, variable names, and actual coding follow those 
used in the programs for potential problems, particularly the quadratic element program described in 
Chapter 4. As in the case of that program, the principle adopted in programming is to try to make the 
coding straightforward to understand and follow, rather than necessarily the most efficient in terms 
of computation. 
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6.1 Program BEM2EQ 

The program name indicates that it is for boundary element method analysis of two-dimensional elastic 
problems using quadratic elements. Each of the subprogram units which make up the whole are described 
in turn. The Preface explains how the full program can be accessed as a single file. 

As with programs for potential problems, much of the coding is concerned with the definition of the 
arrangement of elements on the boundary (or boundaries) of the solution domain, and the application of 
boundary conditions to them. Again, entering the co-ordinates of the nodes of each and every element, 
followed by the type and magnitude of the relevant boundary condition applied to each is an option, 
but tedious. Instead, each boundary can be divided up into a series of boundary segments, which are 
either straight lines or circular arcs. The number of elements within a segment can then be chosen, and 
varied easily, and from which the program generates all the element geometric data. The elements on 
a segment do not have to be uniform in size, but can be varied in length by a constant ratio between 
successive elements. In the present version of the program, each segment is subject to only one uniform 
boundary condition: defined displacements or defined stresses. The program distributes this condition to 
all the elements involved. Consequently, the ends of segments are conveniently defined as points where 
there is a significant change in either shape (a corner, for example) or boundary condition. Use of this 
facility is demonstrated later in this chapter. 

6.1.1 Main program 

At the beginning of the program is a storage module named SHAREDDATA2EQ which allows stored data 
to be accessed and shared by all those subprograms that require them (by means of a USE statement). 
The dimensioned array sizes in the module allow for up to 250 quadratic boundary elements with 500 
nodes, with up to 10 different boundaries forming the solution domain, and a maximum number of 
point (node) displacement constraints of 20. With two equations per node, this means that up to 1000 
equations can be solved. A dictionary of the variable names used is provided at the beginning of the 
book, and at the beginning of Part II. 

The main program named BEM2EQ is designed mainly to call each of the other subprograms in turn. It 
does, however, also serve to name the files with which the program communicates via OPEN statements. 
File DATA, which is addressed as file number 5 in the program, serves to supply the input data which 
defines the current problem. The main output of results is to file RESULTS, addressed as file number 
6. Element mesh data, on the other hand, are output to file MESHRES (mesh results) and numbered 7. 

MODULE SHAREDDATA2EQ 

i 

! MODULE STORING SHARED DATA. 

i 
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REAL :: XEEND(260),YEEND(260),XNODE(500),YNODE(500) 

REAL :: XSEND(250),YSEND(250) 

REAL :: UNMX(250,3),UNMY(250,3),ANGSTORE(260),MAXL 

REAL :: A(1000,1001),UV(1000),U(500),V(500),USEG(250),VSEG(250) 

REAL :: AROWX(250,6),AROWY(250,6),BROWX(250,6),BROWY(250,6) 

REAL :: ZG(8),WG(8),EGL(8),WGL(8),JACOB,UNGX,UNGY 

REAL :: SIGNNSEG(250),SIGSNSEG(250),UELEM(250),VELEM(250) 

REAL :: SIGNN(250,3),SIGSS(250,3),SIGSN(250,3),SIGE(250,3) 

REAL :: TX(500),TY(500),TMX(250,3),TMY(250,3) 

REAL :: FXSEG(250),FYSEG(250) 

REAL :: SF(3,8),SD(3,13),SFL(4,3,8),SDL(4,3,8) 

REAL :: PI,E,NU,ESTORE 

REAL :: AKXX,AKXY,AKYX,AKYY,BKXX,BKXY,BKYX,BKYY 
REAL :: BK2XX,BK2XY,BK2YX,BK2YY 

INTEGER :: NEL,NNP,MAXNEL,MAXNNP,MAXNB,NEEND,NGAUSS 
INTEGER :: NODE(250,3),Ml(250),M3(250) 

INTEGER :: NBOUND,NSEGTOT,NELB(10),NSEGB(10) 

INTEGER :: NBCU,NBCS,NBCM,NBCT,IBCE(250),IBCN(500),ISEGBC(250) 
INTEGER :: ISEGEND(250),ISEGELEM(250),MFIRST(250),MLAST(250) 
INTEGER :: NBCPC,MAXNPC,NODEPC(20),IDIRPC(20) 
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END MODULE SHAREDDATA2EQ 


PROGRAM BEM2EQ 

i 

! PROGRAM FOR SOLVING TWO DIMENSIONAL ELASTIC STRESS ANALYSIS PROBLEMS 
! BY THE BOUNDARY ELEMENT METHOD USING QUADRATIC ELEMENTS. 

i 

USE SHAREDDATA2EQ 
OPEN(5,FILE="DATA") 

OPEN(6,FILE="RESULTS") 

OPEN(7,FILE="MESHRES") 

PI=4.0*ATAN(1.) 

i 

! DEFINE THE MAXIMUM PROBLEM SIZE PERMITTED BY THE ARRAY DIMENSIONS. 
MAXNEL=250 
MAXNNP=500 
MAXNB=10 
MAXNPC=2 0 

i 

! INPUT THE PROBLEM TITLE AND TYPE, ALSO MATERIAL PROPERTIES. 

CALL INTITLE 

i 

! INPUT AND GENERATE THE MESH DATA. 

CALL MESHQ 

i 

! OUTPUT THE MESH DATA. 

CALL MSHOUT 

i 

! EVALUATE AND STORE VALUES OF THE SHAPE FUNCTIONS AND THEIR DERIVATIVES 
! AT THE GAUSS POINTS AND NODES. 

CALL SHAPE 

i 

! INPUT, PROCESS AND OUTPUT THE BOUNDARY CONDITIONS. 

CALL BCS 

i 

! FORM THE COEFFICIENT MATRIX AND APPLY THE BOUNDARY CONDITIONS. 

CALL FRMTRX 
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! SOLVE THE LINEAR EQUATIONS. 

NEQN=2 *NNP 
MAXNEQN=2*MAXNNP 
MAXNEQNP1=MAXNEQN+1 

CALL ELIMIN(A,UV,NEQN,MAXNEQN,MAXNEQNPl,IFLAG) 

IF(IFLAG == 1) THEN 
WRITE(6,61) 

61 FORMAT(/ "MATRIX ILL-CONDITIONING DETECTED IN EQUATION SOLVER") 

STOP 
END IF 

i 

! OUTPUT THE NODAL DISPLACEMENTS, ELEMENT STRESSES AND FORCES ON THE 
! BOUNDARY SEGMENTS. 

CALL OUTPUT 

i 

STOP 

END PROGRAM BEM2EQ 

After defining the maximum numbers of boundary elements (MAXNEL), nodes (MAXNNP), boundaries 
(MAXNB) and point constraints (MAXNPC) permitted by the array dimensions, the main program calls 
the following subprograms in turn: 

INTITLE for the problem title, 

MESHQ to input and create the mesh data, 

MSHOUT to write out the mesh data, 

SHAPE for defining the element shape functions, 

BCS for the boundary conditions, 

FRMTRX to define the [A] and [B] coefficient matrices, 

ELIMIN to solve the equations, and finally 
OUTPUT to write out the results. 

As in the potential programs, the matrix ill-conditioning warning is most likely to be triggered by trying 
to solve a poorly defined problem, such as one with only prescribed stress boundary conditions, and 
displacement nowhere defined. Computed results from ELIMIN are contained in array UV: typically 
nodal point displacements, but tractions for nodes subject to prescribed displacements. 
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6.1.2 Subprogram INTITLE 

An alphanumeric title for the problem is first read into TITLE. Next a six character message, either 
“STRESS” or “STRAIN” is read into CASE, to define whether the domain is in a state of plane stress or 
plane strain. The default is plane stress: any message other than STRAIN will result in a state of plane 
stress being assumed. Finally, the values of Youngs modulus and Poissons ratio are read into E and 
NU. If plane strain is to be assumed, the effective values of these properties are defined according to 
Equations 5.7. 

SUBROUTINE INTITLE 

i 

! SUBPROGRAM TO INPUT PROBLEM TITLE AND WHETHER PLANE 
! STRESS OR PLANE STRAIN. ALSO THE MATERIAL PROPERTIES. 

j 

USE SHAREDDATA2EQ 
CHARACTER(80) :: TITLE 

CHARACTER(6) :: CASE 

j 

! INPUT THE PROBLEM TITLE. 

READ(5,FMT="(A80)") TITLE 
WRITE(6,61) TITLE 
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61 FORMAT("QUADRATIC BOUNDARY ELEMENT SOLUTION FOR", 

& " TWO DIMENSIONAL ELASTIC PROBLEM" // A) 

i 

! INPUT THE PROBLEM CASE TYPE: 

! "STRESS" DEFINES PLANE STRESS 

! "STRAIN" DEFINES PLANE STRAIN 

! THE DEFAULT IS PLANE STRESS. 

READ(5,FMT="(A6)") CASE 
WRITE(6,62) CASE 

62 FORMAT(/"UNDER PLANE ",A," CONDITIONS") 

i 

! INPUT AND OUTPUT YOUNG'S MODULUS AND POISSON'S RATIO. 

READ(5,*) E,NU 
WRITE(6,63) E,NU 

63 FORMAT(/"YOUNG'S MODULUS = ",E12.4,1OX,"POISSON'S RATIO = ",F5.3) 

i 

! MODIFY PROPERTIES IF CASE IS ONE OF PLANE STRAIN. 

IF(CASE == "STRAIN") THEN 
E=E/(1.-NU* *2) 

NU=NU/(1.-NU) 

END IF 

i 

RETURN 

END SUBROUTINE INTITLE 

6.1.3 Subprogram to input and generate the element mesh data 

Subprogram MESHQ is identical to the subprogram described in Section 4.1.3, with the exception that 
the USE statement is for module SHAREDDATA2EQ rather than SHAREDDATA2PQ. 

SUBROUTINE MESHQ 

i 

! SUBPROGRAM TO READ IN AND GENERATE THE GEOMETRIC DATA FOR A MESH OF 
! QUADRATIC ELEMENTS. 

I 

USE SHAREDDATA2EQ 

I 

! INPUT THE NUMBER OF SEPARATE BOUNDARIES. 

READ(5,*) NBOUND 
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! TEST THE NUMBER OF BOUNDARIES. 

IF(NBOUND < 1 .OR. NBOUND > MAXNB) THEN 
WRITE(6,61) NBOUND,MAXNB 

61 FORMAT(/ "NBOUND =",14,2X,"OUTSIDE PERMITTED RANGE 1 TO",14) 

STOP 

END IF 

j 

! FOR EACH BOUNDARY IN TURN INPUT THE NUMBER OF SEGMENTS. 

NEL=0 

IEEND=0 

NSEGTOT=0 

MMAXB=0 

Each boundary in turn: DO IB0UND=1,NBOUND 

NELB(IBOUND)=0 

MMINB=MMAXB+1 

READ(5,*) NSEGB(IBOUND) 

NSEGTOT=NSEGTOT+NSEGB(IBOUND) 

j 

! TEST THE NUMBER OF SEGMENTS. 

IF(NSEGTOT < 1 .OR. NSEGTOT > MAXNEL) THEN 
WRITE(6,62) NSEGTOT,MAXNEL 

62 FORMAT(/ "NSEGTOT =",I6,2X,"OUTSIDE PERMITTED RANGE 1 TO",16) 

STOP 

END IF 

j 

! INPUT THE CARTESIAN GLOBAL COORDINATES OF THE END POINTS OF THE 
! SEGMENTS. TAKE THE END POINTS CONSECUTIVELY, KEEPING THE DOMAIN 
! TO THE LEFT OF THE DIRECTION OF NUMBERING. 

READ(5,*) (XSEND(ISEND),YSEND(ISEND),ISEND=1,NSEGB(IBOUND)) 

j 

! DEFINE THE FIRST END POINT ON THE CURRENT BOUNDARY. 

IEEND=IEEND+1 

XEEND(IEEND)=XSEND(1) 

YEEND(IEEND)=YSEND(1) 

| 

! FOR EACH OF THE SEGMENTS (BETWEEN ENDS 1 AND 2, 2 AND 3, ETC.) 

! INPUT THE RADIUS OF CURVATURE (+VE FOR CONVEX WITH CENTRE OF 
! CURVATURE INSIDE DOMAIN, -VE FOR CONCAVE), THE NUMBER OF 
! ELEMENTS IN THE SEGMENT, AND THE LENGTH RATIO BETWEEN SUCCESSIVE 
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! ELEMENTS IN THE DIRECTION OF NUMBERING. 

ISEGMAX=NSEGTOT 

ISEGMIN=ISEGMAX-NSEGB(IBOUND)+1 

Each segment in turn: DO ISEG=ISEGMIN,ISEGMAX 

READ(5, *) RSEG,NELSEG,RATSEG 

i 

! FIND AND TEST THE NUMBER OF ELEMENTS SO FAR. 

NEL=NEL+NELSEG 

NELB(IBOUND)=NELB(IBOUND)+NELSEG 
IF(NEL <1 .OR. NEL > MAXNEL) THEN 
WRITE(6,63) NEL,MAXNEL 

63 FORMAT(/ "NEL =",16,2X,"OUTSIDE PERMITTED RANGE 1 TO",16) 

STOP 
END IF 

i 

! FIRST AND LAST ELEMENTS ON CURRENT SEGMENT. 

MLAST(ISEG)=NEL 

MFIRST(ISEG)=NEL-NELSEG+1 

MMAXB=MMAXB+NE L S E G 
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! COORDINATES OF THE FIRST END POINT OF THE SEGMENT. 
ISEND=ISEG-ISEGMIN+1 
XFIRST=XSEND(ISEND) 

YFIRST=YSEND(ISEND) 

j 

! COORDINATES OF THE LAST END POINT OF THE SEGMENT. 

ISEND=ISEND+1 

IF(ISEG == ISEGMAX) ISEND=1 
XLAST=XSEND(ISEND) 

YLAST=YSEND(ISEND) 

| 

! GENERATE ELEMENT DATA FOR A STRAIGHT SEGMENT. 

IF(RSEG == 0.) THEN 

j 

! DEFINE THE ELEMENT END POINT COORDINATES ON THE SEGMENT. 

Each element in turn: DO M=1,NELSEG 
IEEND=IEEND+1 
ISEGEND(IEEND)=ISEG 
IF(RATSEG == 1.) THEN 

XEEND(IEEND)=XFIRST+(XLAST-XFIRST)*FLOAT(M)/FLOAT(NELSEG) 
YEEND(IEEND)=YFIRST+(YLAST-YFIRST)*FLOAT(M)/FLOAT(NELSEG) 
ENDIF 

IF(RATSEG /= 1.) THEN 

XEEND(IEEND)=XFIRST+(XLAST-XFIRST)*(1.-RATSEG**M) 

& /(1.-RATSEG**NELSEG) 

YEEND(IEEND)=YFIRST+(YLAST-YFIRST)*(1.-RATSEG**M) 

& /(1.-RATSEG**NELSEG) 

END IF 

END DO Each element in turn 

j 

! DEFINE THE ELEMENT NODES AND COORDINATES. 

IEEND=IEEND-NELSEG 

Each element in turn: DO IELSEG=1,NELSEG 

IEEND=IEEND+1 

M=MFIRST(ISEG)+IELSEG-1 

I1=2*M-1 

12 = 11+1 

13=11+2 

IF(ISEG == ISEGMAX .AND. IELSEG == NELSEG) I3=NODE(MMINB,1) 
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NODE(M,1)=11 
NODE(M,2)=12 
NODE(M,3)=13 
ISEGELEM(M)=ISEG 

IF(ISEG == ISEGMIN .AND. IELSEG == 1) THEN 
XNODE(II)=XEEND(IEEND-1) 

YNODE(II)=YEEND(IEEND-1) 

END IF 

XNODE(13)=XEEND(IEEND) 

YNODE(13)=YEEND(IEEND) 

XNODE(12)=0.5*(XNODE(II)+XNODE(13)) 

YNODE(12)=0.5*(YNODE(II)+YNODE (13) ) 

j 

! STORE THE NUMBERS OF THE ADJACENT ELEMENTS. 

Ml(M)=M-1 
M3(M)=M+1 

END DO Each element in turn 
END IF 

j 

! GENERATE ELEMENT DATA FOR A SEGMENT IN THE FORM OF A CIRCULAR ARC. 

IF(RSEG /= 0.) THEN 

| 

! LOCATE THE CENTRE OF THE ARC. 

XMID=(XFIRST+XLAST)/2. 

YMID=(YFIRST+YLAST)/2. 

ALSEG=SQRT((XLAST-XFIRST)**2+(YLAST-YFIRST)**2) 

ALPERP2=RSEG* *2-(ALSEG/2.) **2 

IF(ABS(ALPERP2) < 1.E-6*RSEG**2) ALPERP2=0. 

IF(ALPERP2 < -1.E-6*RSEG**2) THEN 
WRITE(6,64) ISEG 

64 FORMAT(/ "DATA ERROR FOR SEGMENT NUMBER",16, 

& / "NOT POSSIBLE TO CREATE A CIRCULAR ARC") 

STOP 
END IF 

ALPERP=SQRT(ALPERP2) 

UVFLX=(XLAST-XFIRST)/ALSEG 
UVFLY=(YLAST-YFIRST)/ALSEG 
FACT=1. 

IF(RSEG < 0.) FACT=-1. 
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XCENT=XMID-ALPERP*UVFLY*FACT 

YCENT=YMID+ALPERP*UVFLX*FACT 

| 

! FIND THE ANGLE SUBTENDED THERE BY THE SEGMENT. 

IF(ALPERP /= 0.) ANGSEG=2.*ATAN(ALSEG*0.5/ALPERP) 

IF(ALPERP == 0.) ANGSEG=PI 

j 

! DEFINE THE ELEMENT END POINT COORDINATES ON THE SEGMENT. 

ANGFIR=ATAN2(YFIRST-YCENT,XFIRST-XCENT) 

ANGSTORE(IEEND)=0. 

Each element in turn: DO M=1,NELSEG 

IEEND=IEEND+1 

ISEGEND(IEEND)=ISEG 

IF(RATSEG mm. 1.) ANG=ANGSEG*FLOAT(M)/FLOAT(NELSEG) 

IF(RATSEG /= 1.) ANG=ANGSEG*(1.-RATSEG**M)/(1.-RATSEG**NELSEG) 

IF(RSEG < 0.) ANG=-ANG 

XEEND(IEEND)=XCENT+ABS(RSEG)*COS(ANGFIR+ANG) 

YEEND(IEEND)=YCENT+ABS(RSEG)*SIN(ANGFIR+ANG) 

ANGSTORE(IEEND)=ANG 

END DO Each element in turn 
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i 

! DEFINE THE ELEMENT NODES AND COORDINATES. 

IEEND=IEEND-NELSEG 

Each element in turn: DO IELSEG=1,NELSEG 

IEEND=IEEND+1 

M=MFIRST(ISEG)+IELSEG-1 

11=2 *M-1 

12 = 11+1 

13=11+2 

IF(ISEG == ISEGMAX .AND. IELSEG == NELSEG) I3=N0DE(MMINB,1) 

NODE(M,1)=11 
NODE(M,2)=12 
NODE(M,3)=13 
ISEGELEM(M)=ISEG 

IF(ISEG == ISEGMIN .AND. IELSEG == 1) THEN 
XNODE(II)=XEEND(IEEND-1) 

YNODE(II)=YEEND(IEEND-1) 

END IF 

XNODE(13)=XEEND(IEEND) 

YNODE(13)=YEEND(IEEND) 

ANG=0.5*(ANGSTORE(IEEND-1)+ANGSTORE(IEEND)) 

XNODE(12)=XCENT+ABS(RSEG)*COS(ANGFIR+ANG) 

YNODE(12)=YCENT+ABS(RSEG)*SIN(ANGFIR+ANG) 

I 

! STORE THE NUMBERS OF THE ADJACENT ELEMENTS. 

Ml(M)=M-1 
M3(M)=M+1 

END DO Each element in turn 
END IF 

I 

END DO Each segment in turn 

I 

! ADJACENT ELEMENTS FOR END ELEMENTS OF THE BOUNDARY. 

Ml(MMINB)=MMAXB 

M3(MMAXB)=MMINB 

END DO Each boundary in turn 

NEEND=IEEND 

NNP=NEL*2 
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! DETERMINE THE MAXIMUM DIMENSION OF THE SOLUTION DOMAIN. 

MAXL=0. 

Each node in turn: DO 1=1,NNP 
Each other node in turn: DO J=1,NNP 

DIST=SQRT((XNODE(I)-XNODE(J))**2+(YNODE(I)-YNODE(J))**2) 

IF(DIST > MAXL) MAXL=DIST 
END DO Each other node in turn 
END DO Each node in turn 

I 

RETURN 

END SUBROUTINE MESHQ 

6.1.4 Mesh data output subprogram 

Subprogram MSHOUT is identical to the subprogram described in Section 4.1.4, with the exception that 
the USE statement is for module SHAREDDATA2EQ rather than SHAREDDATA2PQ. 

SUBROUTINE MSHOUT 

j 

! SUBPROGRAM TO WRITE OUT THE MESH DATA. 

| 

USE SHAREDDATA2EQ 

j 

! OUTPUT THE NUMBERS OF ELEMENTS AND NODES, ALSO THE NODAL 
! COORDINATES. 

WRITE(7,71) NEL,NNP,(I,XNODE(I),YNODE(I),1=1,NNP) 

71 FORMAT(/ "GEOMETRIC DATA FOR THE MESH" // 

& 10X,"NUMBER OF ELEMENTS =",I6 // 

& 10X,"NUMBER OF NODAL POINTS =",I6 // 

& "COORDINATES OF THE NODAL POINTS"// 

& 2 (" I X Y ") / 

& 2 (16,2E12.4)) 

| 

! OUTPUT THE ELEMENT NODE NUMBERS. 

WRITE(7,72) (M,(NODE(M,IN),IN=1,3),M=1,NEL) 

72 FORMAT(/ 'ELEMENT NODE NUMBERS' // 

& IX,2(10X,'M ND1 ND2 ND3')/ (2(7X,4I5))) 

| 

! SCALE THE NODAL POINT COORDINATES. 

Each node in turn: DO 1=1,NNP 
XNODE(I)=XNODE(I)/MAXL 
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YNODE(I)=YNODE(I)/MAXL 
END DO Each node in turn 

i 

RETURN 

END SUBROUTINE MSHOUT 

6.1.5 Subprogram for defining shape functions 

Subprogram SHAPE is identical to the subprogram described in Section 4.1.6, with the exception that 
the USE statement is for module SHAREDDATA2EQ rather than SHAREDDATA2PQ. 


SUBROUTINE SHAPE 

I 

! SUBPROGRAM TO EVALUATE AND STORE VALUES OF THE SHAPE FUNCTIONS AND 
! THEIR DERIVATIVES, AT THE GAUSS POINTS AND NODES. 

i 

USE SHAREDDATA2EQ 

i 

! STORE APPROPRIATE COORDINATES AND WEIGHT FACTORS FOR NORMAL GAUSSIAN 
! QUADRATURE IN ARRAYS XG AND CG, ALSO THOSE FOR LOGARITHMIC FUNCTION 
! INTEGRATION IN XGL AND CGL. 
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NGAUSS=8 

ZG(l)=-0.9602898564 
ZG(2)=-0.7966664774 
ZG(3)=-0.5255324099 
ZG(4)=-0.1834346424 
ZG(5)=-ZG(4) 

ZG(6)=-ZG(3) 

ZG(7)=-ZG(2) 

ZG(8)=-ZG(1) 

WG(1)=0.1012285362 
WG(2)=0.2223810344 
WG(3)=0.3137066458 
WG(4)=0.3626837833 
WG(5)=WG(4) 

WG(6)=WG(3) 

WG(7)=WG(2 ) 

WG(8)=WG(1) 

EGL(1)=0.013320244 
EGL(2)=0.079750429 
EGL(3)=0.197871029 
EGL(4)=0.354153994 
EGL(5)=0.529458575 
EGL(6)=0.701814530 
EGL(7)=0.849379320 
EGL(8)=0.953326450 
WGL(1)=0.164416605 
WGL(2)=0.237525610 
WGL(3)=0.226841984 
WGL(4)=0.175754079 
WGL(5)=0.112924030 
WGL(6)=0.057872211 
WGL(7)=0.020979074 
WGL(8)=0.003686407 

i 

! NORMAL GAUSSIAN QUADRATURE. 

For each Gauss point in turn: DO IGAUSS=1,NGAUSS 
ZETA=ZG(IGAUSS) 

SF(1 , IGAUSS)=0.5*ZETA*(ZETA-1.) 

SF(2,IGAUSS) =1 .-ZETA* *2 

SF(3,IGAUSS)=0.5*ZETA*(ZETA+1.) 
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SD(1,IGAUSS)=ZETA-0.5 
SD(2,IGAUSS)=-2.*ZETA 
SD(3, IGAUSS)=ZETA+0.5 
END DO For each Gauss point in turn 

j 

! FOUR CASES OF LOGARITHMIC GAUSSIAN QUADRATURE TO CONSIDER. 
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TO 
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For each case in turn: DO IC=1,4 

For each Gauss point in turn: DO IGAUSS=1,NGAUSS 
ETA=EGL(IGAUSS) 

IF(IC == 1) ZETA=2.*ETA-1. 

IF(IC == 2) ZETA=ETA 

IF(IC == 3) ZETA=-ETA 

IF(IC == 4) ZETA=1.-2.*ETA 

SFL(IC,1,IGAUSS)=0.5*ZETA*(ZETA-1.) 

SFL(IC,2,IGAUSS)=1.-ZETA**2 

SFL(IC, 3,IGAUSS)=0.5*ZETA*(ZETA+1.) 

SDL(IC,1,IGAUSS)=ZETA-0.5 

SDL(IC,2,IGAUSS)=-2.*ZETA 

SDL(IC,3,IGAUSS)=ZETA+0.5 

END DO For each Gauss point in turn 

END DO For each case in turn 

I 

! SHAPE FUNCTION DERIVATIVES AT THE NODES, STORED AS THOUGH THEY 
! ARE FOR GAUSS POINTS NUMBERED 11, 12 AND 13. 

Each element node in turn: DO IGAUSS=11,13 
IF (IGAUSS == 11) ZETA—1. 

IF(IGAUSS == 12) ZETA=0. 

IF(IGAUSS == 13) ZETA=1. 

SD(1,IGAUSS)=ZETA-0.5 
SD(2,IGAUSS)=-2.*ZETA 
SD(3,IGAUSS)=ZETA+0.5 
END DO Each element node in turn 

j 

RETURN 

END SUBROUTINE SHAPE 
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6.1.6 Subprogram for applying the boundary conditions 

The subprogram BCS serves to apply boundary conditions of either the prescribed displacement or 
prescribed stress types. As already indicated, it is assumed that each segment of elements has a uniform 
boundary condition applied to it, which is an important consideration when defining the segments. Also 
applied are any point displacement constraints. The total numbers of segments subject to the first two 
types of condition, and the number of point constraints, are first read into variables NBCU, NBCS and 
NBCPC, respectively. In the case of prescribed stresses it is only necessary to include those segments 
subject to non-zero stresses. 

Prior to reading in the boundary conditions, the components in the global co-ordinate directions of the 
unit normals at every node are calculated with the aid of subprogram JACOBI. Because the direction of 
the normal to the boundary may be discontinuous at an element end node, particularly at a corner of 
the domain, the calculation is performed for every node of every element, and stored in arrays UNMX 
and UNMY. 
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Storage arrays for normal and shear stresses, SIGNN and SIGSN, and tractions in the global co-ordinate 
directions, TMY and TMY, for every node of every element are set to zero. Array IBCE stores the type 
of boundary condition (1 or 2 for the two types, prescribed displacement or prescribed stress) for each 
element. In the absence of other information, the condition at a node is assumed to be zero stresses. The 
default boundary condition type is therefore set as 2, with the corresponding zero values of stresses for 
the elements being stored in arrays TMX and TMY. 

SUBROUTINE BCS 

i 

! SUBPROGRAM TO INPUT, PROCESS AND OUTPUT THE BOUNDARY CONDITIONS. 

i 

USE SHAREDDATA2EQ 

i 

! FIRST FIND THE UNIT NORMAL COMPONENTS AT THE ELEMENT NODES. 

Each element in turn: DO M=1,NEL 

Each element node in turn: DO IN=1,3 

IGAUSS=IN+10 

IT=1 

IC=1 

CALL JACOBI(M,IGAUSS,IT,IC) 

UNMX(M,IN)=UNGX 
UNMY(M,IN)=UNGY 

END DO Each element node in turn 
END DO Each element in turn 

I 

! INPUT THE NUMBERS OF SEGMENTS SUBJECT TO EACH TYPE OF BOUNDARY 
! CONDITION. ALSO THE NUMBER OF POINT CONSTRAINTS. 

! NBCU - PRESCRIBED DISPLACEMENTS. 

! NBCS - NON-ZERO PRESCRIBED STRESSES. 

! ANY SEGMENT NOT INCLUDED IS ASSUMED TO BE SUBJECT TO ZERO 
! STRESSES. 

! NBCPC - NODAL POINT DISPLACEMENT CONSTRAINTS. 

i 

READ(5,*) NBCU,NBCS,NBCPC 

i 

! TEST THESE BOUNDARY CONDITION NUMBERS. 

NBCT=NBCU+NBCS 

IF(NBCU <0 .OR. NBCU > MAXNEL .OR. NBCS <0 .OR. NBCS > MAXNEL 
& .OR. NBCT <0 .OR. NBCT > MAXNEL) THEN 

WRITE(6,61) NBCU,NBCS,NBCT,MAXNEL 
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61 FORMAT(/ "NBCU = ",16,3X,"NBCS =",I6,3X,/ "NBCT =",I6,3X, 

& "OUTSIDE PERMITTED RANGE 0 TO",16) 

STOP 
END IF 

IF(NBCPC <0 .OR. NBCPC > MAXNPC) THEN 
WRITE(6,62) NBCPC,MAXNPC 

62 FORMAT(/ "NBCPC =",I4,3X,"OUTSIDE PERMITTED RANGE 0 TO",13) 

STOP 

END IF 

i 

! INITIALISE THE BOUNDARY CONDITION STORAGE ARRAYS. 

Each element in turn: DO M=1,NEL 
IBCE(M)=2 

Each element node in turn: DO IN=1,3 
SIGNN(M,IN)=0. 

SIGSN(M,IN)=0. 

TMX(M,IN)=0. 

TMY(M,IN)=0. 

END DO Each element node in turn 
END DO Each element in turn 

i 

! INPUT, STORE AND OUTPUT THE PRESCRIBED DISPLACEMENT CONDITIONS. 

IF(NBCU > 0) THEN 

READ(5,*) (ISEGBC(IBCU),USEG(IBCU),VSEG(IBCU),IBCU=1,NBCU) 

WRITE(6,63) 

63 FORMAT(/ "PRESCRIBED DISPLACEMENT BOUNDARY CONDITIONS") 

Each segment with prescribed displacements: DO IBCU=1,NBCU 
ISEG=ISEGBC(IBCU) 

IF(ISEG < 1 .OR. ISEG > NSEGTOT) THEN 
WRITE(6,64) ISEG,NSEGTOT 

64 FORMAT(/ "ISEG = ",16,2X,"OUTSIDE PERMITTED RANGE 1 TO",16) 

STOP 

END IF 

Each element on current segment: DO M=MFIRST(ISEG),MLAST(ISEG) 
IBCE(M)= 1 

UELEM(M)=USEG(IBCU) 

VELEM(M)=VSEG(IBCU) 

END DO Each element on current segment 

i 

WRITE(6,65) USEG(IBCU),VSEG(IBCU),MFIRST(ISEG),MLAST(ISEG) 
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65 FORMAT(/ "U =",E12.4, 3X, "V = ",E12.4,3X,"ON ELEMENTS",14, 3X, 

& "TO",14) 

END DO Each segment with prescribed displacements 
END IF 

i 

! INPUT, STORE AND OUTPUT THE PRESCRIBED STRESS CONDITIONS. 

IF(NBCS > 0) THEN 

READ(5,*) (ISEGBC(IBCS),SIGNNSEG(IBCS),SIGSNSEG(IBCS), 

& IBCS=1,NBCS) 

WRITE(6,66) 

66 FORMAT(/ "PRESCRIBED STRESS BOUNDARY CONDITIONS") 

Each segment with prescribed stresses: DO IBCS=1,NBCS 
ISEG=ISEGBC(IBCS) 

IF(ISEG < 1 .OR. ISEG > NSEGTOT) THEN 
WRITE(6,64) ISEG,NSEGTOT 
STOP 
END IF 

Each element on current segment: DO M=MFIRST(ISEG),MLAST(ISEG) 
IBCE(M)=2 
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! FIND THE TRACTIONS AT THE ELEMENT NODES FROM THE PRESCRIBED STRESSES. 

! ALSO STORE THE STRESSES AT THE ELEMENT NODES. 

Each element node in turn: DO IN=1,3 

TMX(M,IN)=SIGNNSEG(IBCS)*UNMX(M,IN)-SIGSNSEG(IBCS)*UNMY(M,IN) 

TMY(M,IN)=SIGNNSEG(IBCS)*UNMY(M,IN)+SIGSNSEG(IBCS)*UNMX(M,IN) 
SIGNN(M,IN)=SIGNNSEG(IBCS) 

SIGSN(M,IN)=SIGSNSEG(IBCS) 

END DO each element node in turn 

END DO Each element on current segment 

i 

WRITE(6,67) SIGNNSEG(IBCS),SIGSNSEG(IBCS),MFIRST(ISEG), 

& MLAST(ISEG) 

67 FORMAT(/ "SIGNN =",E12.4,3X,"SIGSN =",E12.4,3X,"ON ELEMENTS", 

& I4, 3X,"TO",14) 

END DO Each segment with prescribed stresses 
END IF 

i 

! ASSEMBLE THE VECTOR OF KNOWN VARIABLES, APPLYING SCALING. 

! FIRST INITIALISE THE NODAL DISPLACEMENTS AND TRACTIONS. 

Each node in turn: DO 1=1,NNP 
U(I)=0. 

V(I)=0. 

TX(I)=0. 

TY(I)=0. 

END DO Each node in turn 
Each element in turn: DO M=1,NEL 
Each element node in turn: DO IN=1,3 
I=NODE(M,IN) 

IBCN(I)=IBCE(M) 

I 

! CHECK THE CONDITION APPLIED TO THE ADJACENT ELEMENT FOR AN ELEMENT 
! END NODE. IF THE CONDITION IS PRESCRIBED DISPLACEMENTS BUT THE 
! EXISTING CONDITION AT THE NODE IS NOT, THEN IMPOSE PRESCRIBED 
! DISPLACEMENTS. 

MADJ=M 

IF(IN == 1) MADJ=M1(M) 

IF(IN == 3) MADJ=M3(M) 

IF(IBCE(MADJ) == 1 .AND. IBCN(I) == 2) IBCN(I)=1 


74 


Download free eBooks at bookboon.com 



Boundary Element Methods for Engineers: 

Part II: Plane Elastic Problems QuadraticBoundaryElementProgramforPlaneElasticProblems 

! STORE KNOWN VARIABLES. 

IF(IBCN(I) == 1) THEN 

IF(IBCE(M) == 1) THEN 

IF(U(I) == 0.) U(I)=UELEM(M)/MAXL 
IF(V(I) == 0.) V(I)=VELEM(M)/MAXL 
IF(U(I) /= 0.) U(I)=0.5*(U(I)+UELEM(M)/MAXL) 

IF(V(I) /= 0.) V(I)=0.5*(V(I)+VELEM(M)/MAXL) 

END IF 
END IF 

IF(IBCE(M) == 2 .AND. IBCN(I) == 2) THEN 
IF(TX(I) == 0.) TX(I)=TMX(M,IN)/E 
IF(TY(I) == 0.) TY(I)=TMY(M,IN)/E 
IF (TX (I) /= 0.) TX (I) =0.5* (TX (I) +TMX (M, IN) /E) 

IF(TY (I) /= 0.) TY (I) =0.5* (TY (I)+TMY (M, IN)/E) 

END IF 

END DO Each element node in turn 
END DO Each element in turn 

j 

! SCALE YOUNG'S MODULUS. 

ESTORE=E 
E=1. 

| 

! INPUT, STORE AND OUTPUT NODAL POINT DISPLACEMENT CONSTRAINTS. 

! IDIRPC STORES DIRECTIONS OF THE CONSTRAINTS - 1 FOR X, 2 FOR Y. 

! SUCH CONSTRAINTS SHOULD NOT REQUIRE POINT FORCES TO MAINTAIN THEM. 

IF(NBCU == 0 .AND. NBCPC < 3) THEN 
WRITE(6,68) 

68 FORMAT(/ "INSUFFICIENT DISPLACEMENT BOUNDARY CONDITIONS", 

& " OR POINT CONSTRAINTS") 

STOP 
END IF 

IF(NBCPC > 0) THEN 

READ(5,*) (NODEPC(IBCPC),IDIRPC(IBCPC),IBCPC=1,NBCPC) 

Each point constraint in turn: DO IBCPC=1,NBCPC 
IF(NODEPC(IBCPC) < 0 .OR. NODEPC(IBCPC) > NNP) THEN 
WRITE(6,69) NODEPC(IBCPC),NNP 

69 FORMAT(/ "POINT CONSTRAINT NODE",14," OUTSIDE RANGE 0 TO ",I4) 
STOP 

END IF 

IF(IDIRPC(IBCPC) / = 1 .AND. IDIRPC(IBCPC) / = 2) THEN 
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70 


& 


WRITE(6,70) IDIRPC(IBCPC),NODEPC(IBCPC) 

FORMAT(/ "POINT CONSTRAINT DIRECTION",I3," FOR NODE",14, 
" OUTSIDE RANGE 1 TO 2") 


71 

72 


STOP 
END IF 

IF(IDIRPC(IBCPC) = 
FORMAT(/ "NODE",14 
IF(IDIRPC(IBCPC) = 
FORMAT(/ "NODE",14 


1) WRITE(6,71) NODEPC 
CONSTRAINED WITH U = 

2) WRITE(6,72) NODEPC 
CONSTRAINED WITH V = 


(IBCPC) 

0 ") 

(IBCPC) 

0 ") 


END DO Each point constraint in turn 


END IF 


RETURN 

END SUBROUTINE BCS 


For each segment over which displacements are prescribed, the segment number and values of 
displacements are read into ISEGBC, USEG and VSEG, respectively. Taking each of these segments 
in turn, the elements have their boundary condition type set to 1 (in array IBCE), and the values of 
displacements are stored in arrays UELEM and VELEM. The prescribed displacements and numbers of 
the elements to which they are applied are then written out. 


For each segment over which stresses are prescribed, the segment number and values of normal and shear 
stresses are read into ISEGBC, SIGNNSEG and SIGSNSEG, respectively. Taking each of these segments 
in turn, the elements have their boundary condition type set to 2 (in array IBCE). The values of the 
stresses are both stored as nodal point stresses in arrays SIGNN and SIGNS and converted to tractions 
at the nodes of the elements according to Equations 5.58 and 5.59 and stored in arrays TMX and TMY. 
The prescribed stresses and numbers of the elements to which they are applied are then written out. 


For each node in turn, the known variable values are stored: displacements in arrays U and V, or tractions 
in arrays TX and TY. Displacements are scaled by dividing by MAXL, the maximum domain dimension. 
Tractions are scaled by dividing by Youngs modulus. The type of boundary condition at the node is 
stored in array IBCN, obtained from that for the corresponding element. At the ends of segments, which 
may correspond to physical corners in the problem or may correspond only to changes in boundary 
conditions, this can give rise to an ambiguity in type at a node. At a node at which there is a change from 
displacement to stress boundary condition the former is chosen, which is why in the program there is a 
test for displacement boundary condition on the current element, but stress boundary condition at the 
current node (which can arise from the way the boundary conditions have been applied to the elements). 
The displacement boundary condition type is imposed by setting the relevant value in array IBCN to 1. 
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If at a node linking two segments there is a change, not in boundary condition type but in the values of 
either displacements or stresses, it is appropriate to assign average values to the node and store these as 
the known variables for that node. This is achieved in the program by testing every node to see whether 
both the node and element boundary condition is of the same type. If it is, the known variables are 
taken as the stored values for the element if the known variable values are zero. This will apply to all 
element midside nodes, and to all element end nodes which have not been tested in previous elements. 
If the known variable values are not zero, which can only occur if prescribed values are applied to the 
immediately adjacent element and that conditions at the nodes of the element have already been tested, 
then the average of the present element values and the stored values are re-stored. This works both for 
adjacent elements with different values of displacement or tractions and for the much more common 
case of adjacent elements with the same value prescribed: the common value is re-stored. 

This averaging process is really a detail to give sensible values of the variables in the output of the final 
results. What is really important is that the prescribed displacements and tractions are applied element 
by element in forming the [ b ] column vector in Equations 5.57, so that each of the adjacent elements 
contributes according to its own prescribed values. This will be explained further in connection with 
subprogram FRMTRX. 


In preparation for using scaled values of the displacement kernels, the value of Youngs modulus is set 
to 1, having stored the actual value in the variable ESTORE. 
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Finally in subprogram BCS the nodal point displacement constraints are read in: the node number into 
NODEPC and the direction of the constraint into IDIRPC. Direction 1 is the % direction and direction 
2 is the y direction. A value of 1 means that the node is constrained not to move in the x direction, 
2 that it cannot move in the y direction. In the absence of any displacement boundary conditions, at 
least 3 (and often only 3) point constraints are required to prevent movement of the domain as a rigid 
body, without at the same time unnaturally restricting deformation of the elastic component under 
consideration. Typically, one node is constrained in both co-ordinate directions, while another which is 
suitably located is constrained to move only in the direction which is radial to the fully constrained node. 
Before the point constraint data are read in, a test is applied to the numbers of displacement boundary 
conditions and point constraints in the subprogram, followed by a warning message where necessary. 

6.1.7 Subprogram to form the coefficient matrix 

Subprogram FRMTRX forms the coefficient matrix and right hand side vector of Equations 5.57, and 
stores them as an extended matrix in array A. The extended matrix is simply [A*] with an extra column 
to contain column vector [b]. Each node in turn is treated as point P, the force point for the fundamental 
solution, then integration is carried out over every element in turn. This integration involves the knowns 
and unknowns at each of the three nodes of each element. Node counters I and J are used for P and 
the current element node, respectively, so that 2*1-1 and 2*1 are the numbers of the two rows in [A*] 
corresponding to P, while 2* J-l and 2*J are the numbers of the two columns in [A*] corresponding to 
the element node. The actual integrations of the kernel function products are carried out in subprogram 
INTKER which is called by FRMTRX. What INTKER does is to compute the integrals of the traction 
and displacement kernel functions required in Equations 5.78 and 5.79. If the knowns and unknowns 
at the th node of the current element are tractions and displacements, then these integral quantities will 
contribute to the [A] and [B] matrices, respectively. Initially they are stored (in subprogram INTKER) 
in the two-dimensional arrays AROWX and AROWY as the coefficients that will form the current pair 
of rows of [A], and BROWX and BROWY as the coefficients that will form the current pair of rows of 
[5], but with the contributions from every node of every element kept distinct. The two rows forming 
the pair for each node are for concentrated force applied at the force point P in the x direction (AROWX 
and BROWX) and y direction (AROWY and BROWY). 

SUBROUTINE FRMTRX 

i 

! SUBPROGRAM TO FORM THE COEFFICIENT MATRIX AND RIGHT HAND SIDE VECTOR, 

! MODIFIED TO SUIT THE BOUNDARY CONDITIONS. 

i 

USE SHAREDDATA2EQ 

i 

! DEFINE THE NUMBER OF COLUMNS IN THE EXTENDED COEFFICIENT MATRIX A. 

JMAX=2 *NNP+1 

i 
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! FORM THE COEFFICIENT MATRIX A, AND THE RIGHT HAND SIDE VECTOR B*T. 

Take each node in turn as P: DO 1=1,NNP 

I 

! INITIALISE THE TWO ROWS OF MATRIX A CORRESPONDING TO THE NODE. 

Each pair of coefficients in turn: DO J=1,JMAX 
A(2*1-1,J)=0. 

A(2*1,J)=0. 

END DO Each pair of coefficients in turn 

I 

! INITIALISE THE ELEMENT CONTRIBUTIONS TO THE CURRENT ROWS OF THE 
! A AND B MATRICES. 

Each element in turn: DO M=1,NEL 
Each element node in turn: DO IN=1,3 
AROWX(M,2*IN-1)=0. 

AROWX(M,2*IN)=0. 

AROWY(M,2*IN-1)=0. 

AROWY(M,2*IN)=0. 

BROWX(M,2*IN-1)=0. 

BROWX(M,2*IN)=0. 

BROWY(M,2*IN-1)=0. 

BROWY(M,2*IN)=0. 

END DO Each element node in turn 
END DO Each element in turn 

j 

! SET UP THE LOOP TO INTEGRATE OVER EACH ELEMENT IN TURN. 

Each element in turn: DO M=1,NEL 

| 

! INTEGRATE THE KERNEL PRODUCTS OVER THE CURRENT ELEMENT. 

CALL INTKER(I,M) 

END DO Each element in turn 

j 

! EVALUATE THE COEFFICIENTS AT THE DIAGONAL OF MATRIX A. 

AIIXX=0. 

AIIXY=0. 

AIIYX=0. 

AIIYY=0. 

Each element in turn: DO M=1,NEL 
Each element node in turn: DO IN=1,3 
AIIXX=AIIXX-AROWX(M,2*IN-1) 
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AIIXY=AIIXY-AROWX(M,2*IN) 

AIIYX=AIIYX-AROWY(M,2*IN-1) 

AIIYY=AIIYY-AROWY(M,2*IN) 

END DO Each element node in turn 
END DO Each element in turn 
IF(IBCN(I) == 2) THEN 

A(2*1-1,2*1-1)=AIIXX 
A(2*1-1,2*1)=AIIXY 
A(2*I,2*1-1)=AIIYX 
A(2*I,2*1)=AIIYY 
END IF 

j 

! INITIALISE THE B*T VECTOR COEFFICIENTS. 
BTX=0. 

BTY=0. 

IF(IBCN(I) == 1) THEN 

BTX=-AIIXX*U(I)-AIIXY*V(I) 
BTY=-AIIYX*U(I)-AIIYY*V(I) 

END IF 
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! APPLY THE BOUNDARY CONDITIONS TO THE CURRENT ROWS OF A AND B, BY 
! CONSIDERING EACH ELEMENT IN TURN. 

Each element in turn: DO M=1,NEL 
Each element node in turn: DO IN=1,3 
J=NODE(M,IN) 

j 

! IF DISPLACEMENTS ARE PRESCRIBED OVER THE ELEMENT, 

! INTERCHANGE THE A AND B COEFFICIENTS. 

IF(IBCE(M) == 1) THEN 

A(2*I-1,2*J-1)=A(2*1-1,2*J-l)-BROWX(M,2*IN-1) 

A(2*1-1,2*J)=A(2*1-1,2*J)-BROWX(M,2*IN) 

A(2*I,2*J-l)=A(2*I,2*J-l)-BROWY(M,2*IN-1) 

A(2 * 1,2 * J)=A(2 *1,2 * J)-BROWY(M,2*IN) 

BTX=BTX-AROWX(M,2*IN-1)*U(J)-AROWX(M,2*IN)*V(J) 

BTY=BTY-AROWY(M,2*IN-1)*U(J)-AROWY(M,2*IN)*V(J) 

END IF 

| 

! IF STRESSES ARE PRESCRIBED OVER THE ELEMENT, STORE THE A MATRIX 
! COEFFICIENTS, EXCEPT AT A CORNER NODE WHERE PRESCRIBED DISPLACEMENTS 
! HAVE BEEN IMPOSED. 

IF(IBCE(M) == 2) THEN 

BTX=BTX+(BROWX(M,2*IN-1)*TMX(M,IN)+BROWX(M,2*IN)*TMY(M,IN)) 

& /ESTORE 

BTY=BTY+(BROWY(M,2*IN-1)*TMX(M,IN)+BROWY(M,2*IN)*TMY(M,IN)) 

& /ESTORE 

IF(IBCN(J) == 2) THEN 

A(2*1-1,2*J-l)=A(2*1-1,2*J-l)+AROWX(M,2*IN-1) 

A(2*1-1,2*J)=A(2*1-1,2*J)+AROWX(M,2*IN) 

A(2*1,2*J-l)=A(2*1,2*J-l)+AROWY(M,2*IN-1) 

A(2 * 1,2 * J)=A(2 * 1,2 * J)+AROWY(M,2*IN) 

END IF 

IF(IBCN(J) == 1) THEN 

BTX=BTX-AROWX(M,2*IN-1)*U(J)-AROWX(M,2*IN)*V(J) 

BTY=BTY-AROWY(M,2*IN-1)*U(J)-AROWY(M,2*IN)*V(J) 

END IF 
END IF 

j 

END DO Each element node in turn 
END DO Each element in turn 
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! STORE THE B*T COEFFICIENTS AS EXTENSIONS OF MATRIX A. 

A(2*1-1,JMAX)=BTX 
A (2 * I,JMAX)=BTY 

END DO Take each node in turn as P 

i 

! APPLY NODAL POINT DISPLACEMENT CONSTRAINTS. 

IF(NBCPC > 0) THEN 

Each point constraint in turn: DO IBCPC=1,NBCPC 
IROW=2 *(NODEPC(IBCPC)-1)+IDIRPC(IBCPC) 

Each column in turn: DO J=1,JMAX 
A(IROW,J)=0. 

END DO Each column in turn 
A(IROW,IROW)=1. 

END DO Each point constraint in turn 
END IF 

i 

RETURN 

END SUBROUTINE FRMTRX 

At the beginning of the subprogram the current rows of array A and all the coefficients of the AROW and 
BROW arrays are set to zero in anticipation of the assembly process. Then for each boundary element 
in turn subprogram INTKER is called to carry out the integrations. Before any boundary conditions are 
applied, the coefficients of the 2x2 submatrix at the diagonal of the [A] matrix, which include the free 
terms, are computed as explained in connection with Equation 5.80. For each of the four coefficients, the 
sum of corresponding values along the current row of matrix [A] must be zero, which allows the value 
at the diagonal to be equated to the sum of the other values with its sign reversed. 

Off-diagonal coefficients for any element end node have two contributions, from the two elements which 
share it, stored in the relevant locations in AROWX and AROWY. The coefficients on the current pair 
of rows of the product of matrix [B*] with the vector of known variables are to be accumulated in BTX 
and BT Y. These are set to zero initially in anticipation of the element contributions to be added to them. 
If the boundary condition at point P is of the displacement type, the products of the coefficients at the 
diagonal of [A] with the prescribed displacements there are known and must be moved from the left 
hand side to the right hand side of the equation, and stored in BTX and BTY with their signs changed. 

The boundary conditions applied to each of the elements in turn are now used to determine how the terms 
stored in the AROW and BROW arrays must be added to [A] and [B]. It is perhaps helpful to consider 
the relevant fragments of the equations represented by the current rows of the [A] and [ B ] matrices, 
the fragments being for the (2*IN-l) th and the (2*IN) th equations for the element, corresponding to the 
IN th node of the element numbered M, as a result of the integration over the element 
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.... +ARO WX(M,2 * IN-1) x u ; + AROWX(M,2*IN)x Vj .... = 

...,+BROWX(M,2*IN-l)x {t x }j +BROWX(M,2*IN)x {t y } + . ( 6- l) 

.... + AROWY(M,2*IN-l)X Uj+ AROWY(M,2*IN)X Vj .... = 

.... + BROWY(M,2*IN-l)X { t x }j +BROWY(M,2*IN)x [t y } +. (6.2) 

where; is the number of the IN th node of the element (J=NODE(M,IN) in programming terms). These 
parts of a pair of equations are presented in a mixture of program and physical variables in an attempt 
to clarify what is going on. 

Tractions known over the element 

If the tractions over the element are known and the displacements unknown then the above products 
involving BROWX and BROWY, which are known quantities, are added to the right hand side vector 
whose coefficients for the current equation are BTX and BTY. Note that because the traction values 
stored in arrays TMX and TMY have not already been scaled, they are divided by the value of Youngs 
modulus stored in ESTORE. With the displacements at the particular node unknown, on the left hand 
side of the equations the computed AROWX and AROWY values are simply added to the relevant [A] 
matrix coefficients stored in A(I,J). 
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If the displacements at the particular node are known, which implies that the node is an end node at a 
point where the boundary condition changes from displacement type to traction type, then the products 
of the AROWX and AROWY terms with known displacements in Equations 6.1 and 6.2, must be taken 
from the left hand side of the equation to the right hand side and subtracted from BTX and BTY. 

Displacements known over the element 

If the displacements over the element (and hence at all three nodes) are known and the tractions are 
unknown then the terms in Equations 6.1 and 6.2 have to be moved to the opposite sides 

.... -BROWX(M,2*IN-1 )x {t x }j - BROWX(M,2*IN)x [t y } + . = 

.... -AROWX(M,2*IN-l)x Uj - AROWX(M,2*IN)X Vj . (6.3) 

.... - BROWY(M,2*IN-1 )x { t x }j - BROWY(M,2*IN)x { t y } + . = 

.... - AROWY(M,2*IN-l)X Uj- AROWY(M,2*IN)x Vj .... ( 6 . 4 ) 

The products of the AROWX and AROWY terms with the known displacements, which are known 
quantities, are subtracted from the right hand side vector whose coefficients for the current equations 
are BTX and BTY. On the left hand sides of the equations the computed BROWX and BROWY values 
are simply subtracted from the relevant [A] matrix coefficients stored in A(I,J). 

Once assembly of the linear equations for the current point P is complete, the right hand side vector 
coefficients accumulated in BTX and BTY are stored in the last column of extended matrix [A]. The 
assembly process is repeated for all nodes as force point P. 

Finally, any nodal point constraints are applied. As already discussed in Section 5.4, these can only be 
used in positions where they would not give rise to concentrated forces at the points concerned. In other 
words they are used to anchor a domain which is in force equilibrium but would otherwise be free to 
move as rigid body. Point constraints are applied as follows. The relevant equation number to be modified 
is the [2(i — 1) + /c] th , where i is the node number and k is the direction number of the constraint (1 
for x and 2 for y). All the coefficients in that equation, including the one in the right hand side vector, 
are set to zero. The coefficient on the diagonal of matrix [A] is then set to one. These changes have the 
effect of defining the relevant displacement at the node to be zero. 
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6.1.8 Subprogram for integrating kernel function products over an element 

As indicated in the previous section, the purpose of subprogram INTKER is to compute the integrals 
involving the kernel functions for a particular source point P (node numbered I) and element m, and 
which are required in Equations 5.78 and 5.79. The global co-ordinates of P are first stored in XP and 
YP. Then, for each element node in turn (numbered c in the these equations, IN in the subprogram) the 
node number is stored as J. 

SUBROUTINE INTKER(I,M) 

i 

! SUBPROGRAM TO INTEGRATE THE KERNEL PRODUCTS FOR A PARTICULAR FORCE 
! POINT P (INDICATED BY NODE NUMBER I) OVER A PARTICULAR ELEMENT 
! (INDICATED BY ARGUMENT M) BY GAUSSIAN QUADRATURE. 

j 

USE SHAREDDATA2EQ 

| 

! CONSTANT IN DISPLACEMENT KERNELS MULTIPLYING THE LOGARITHMIC TERM. 

CL=(1.+NU)*(3.-NU)/(4.*PI*E) 

j 

! COORDINATES OF POINT P. 

XP=XNODE(I) 

YP=YNODE(I) 

| 

! SET UP THE ELEMENT NODE LOOP. 

Each element node in turn: DO IN=1,3 
J=NODE(M,IN) 

| 

! IF P IS NOT THE CURRENT ELEMENT NODE, USE NORMAL GAUSSIAN QUADRATURE. 
IF(I /= J) THEN 

Each Gauss point in turn: DO IGAUSS=1,NGAUSS 

j 

! EVALUATE JACOBIAN AND UNIT NORMAL VECTOR COMPONENTS, ALSO THE KERNELS 

! AT THE PARTICULAR GAUSS POINT FOR NORMAL QUADRATURE OVER THE WHOLE 
! ELEMENT. 

IT=1 

IC=1 

CALL JACOBI(M,IGAUSS,IT,IC) 

CALL KERNEL(XP,YP,M,IGAUSS,UNGX,UNGY) 

| 

! ACCUMULATE THE INTEGRALS. 

SFN=SF(IN,IGAUSS) 
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AROWX(M,2*IN-1)=AROWX(M,2*IN-1)+WG(IGAUSS)*AKXX*SFN*JACOB 
AROWX(M,2*IN)=AROWX(M,2*IN)+WG(IGAUSS)*AKXY*SFN*JACOB 
AROWY(M,2*IN-1)=AROWY(M,2*IN-1)+WG(IGAUSS)*AKYX*SFN*JACOB 
AROWY(M,2*IN)=AROWY(M,2*IN)+WG(IGAUSS)*AKYY*SFN*JACOB 
BROWX(M,2*IN-1)=BROWX(M,2*IN-1)+WG(IGAUSS)*BKXX*SFN*JACOB 
BROWX(M,2*IN)=BROWX(M,2*IN)+WG(IGAUSS)*BKXY*SFN*JACOB 
BROWY(M,2*IN-1)=BROWY(M,2*IN-1)+WG(IGAUSS)*BKYX*SFN*JACOB 
BROWY(M,2*IN)=BROWY(M,2*IN)+WG(IGAUSS)*BKYY*SFN*JACOB 
END DO Each Gauss point in turn 
END IF 

| 

! IF P IS THE CURRENT ELEMENT NODE, SOME LOGARITHMIC QUADRATURE 
! IS REQUIRED. 

IF(I == J) THEN 

| 

! P AT THE FIRST NODE OF THE ELEMENT. 

IF(IN == 1) THEN 

j 

! TERMS INVOLVING NORMAL QUADRATURE. 

Each Gauss point in turn: DO IGAUSS=1,NGAUSS 
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IT=1 

IC=1 

CALL JACOBI(M,IGAUSS,IT,IC) 

CALL KERN2(M,IN,IGAUSS) 

SFN=SF(IN,IGAUSS) 

BROWX(M,2*IN-1)=BROWX(M,2*IN-1)+WG(IGAUSS)*BK2XX*SFN*JACOB 
BROWX(M,2*IN)=BROWX(M,2*IN)+WG(IGAUSS)*BK2XY*SFN*JACOB 
BROWY(M,2*IN-1)=BROWY(M,2*IN-1)+WG(IGAUSS)*BK2YX*SFN*JACOB 
BROWY(M,2*IN)=BROWY(M,2*IN)+WG(IGAUSS)*BK2YY*SFN*JACOB 
END DO Each Gauss point in turn 

| 

! TERMS INVOLVING LOGARITHMIC QUADRATURE. 

Each Gauss point in turn: DO IGAUSS=1,NGAUSS 

IT=2 

IC=1 

CALL JACOBI(M,IGAUSS,IT,IC) 

SFN=SFL(IC,IN,IGAUSS) 

DZDE=2. 

BROWX(M,2*IN-1)=BROWX(M,2*IN-1)+WGL(IGAUSS)*CL*SFN*JACOB*DZDE 
BROWY(M,2*IN)=BROWY(M,2*IN)+WGL(IGAUSS)*CL*SFN*JACOB*DZDE 
END DO Each Gauss point in turn 
END IF 

| 

! P AT THE SECOND NODE OF THE ELEMENT. 

IF(IN == 2) THEN 

j 

! TERMS INVOLVING NORMAL QUADRATURE. 

Each Gauss point in turn: DO IGAUSS=1,NGAUSS 

IT=1 

IC 1 

CALL JACOBI(M,IGAUSS,IT,IC) 

CALL KERN2(M,IN,IGAUSS) 

SFN=SF(IN,IGAUSS) 

BROWX(M,2*IN-1)=BROWX(M,2*IN-1)+WG(IGAUSS)*BK2XX*SFN*JACOB 
BROWX(M,2*IN)=BROWX(M,2*IN)+WG(IGAUSS)*BK2XY*SFN*JACOB 
BROWY(M,2*IN-1)=BROWY(M,2*IN-1)+WG(IGAUSS)*BK2YX*SFN*JACOB 
BROWY(M,2*IN)=BROWY(M,2*IN)+WG(IGAUSS)*BK2YY*SFN*JACOB 
END DO Each Gauss point in turn 
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! TERMS INVOLVING LOGARITHMIC QUADRATURE. 

Each Gauss point in turn: DO IGAUSS=1,NGAUSS 

IT=2 

IC=2 

CALL JACOBI(M,IGAUSS,IT,IC) 

SFN=SFL(IC,IN,IGAUSS) 

DZDE=1. 

BROWX(M,2*IN-1)=BROWX(M,2*IN-1)+WGL(IGAUSS)*CL*SFN*JACOB*DZDE 
BROWY(M,2*IN)=BROWY(M,2*IN)+WGL(IGAUSS)*CL*SFN*JACOB*DZDE 
IC=3 

CALL JACOBI(M,IGAUSS,IT,IC) 

SFN=SFL(IC,IN,IGAUSS) 

DZDE=1. 

BROWX(M,2*IN-1)=BROWX(M,2*IN-1)+WGL(IGAUSS)*CL*SFN*JACOB*DZDE 
BROWY(M,2*IN)=BROWY(M,2*IN)+WGL(IGAUSS)*CL*SFN*JACOB*DZDE 
END DO Each Gauss point in turn 
END IF 

j 

! P AT THE THIRD NODE OF THE ELEMENT. 

IF(IN == 3) THEN 

| 

! TERMS INVOLVING NORMAL QUADRATURE. 

Each Gauss point in turn: DO IGAUSS=1,NGAUSS 

IT—1 

IC=1 

CALL JACOBI(M,IGAUSS,IT,IC) 

CALL KERN2(M,IN,IGAUSS) 

SFN=SF(IN,IGAUSS) 

BROWX(M,2*IN-1)=BROWX(M,2*IN-1)+WG(IGAUSS)*BK2XX*SFN*JACOB 
BROWX(M,2*IN)=BROWX(M,2*IN)+WG(IGAUSS)*BK2XY*SFN*JACOB 
BROWY(M,2*IN-1)=BROWY(M,2*IN-1)+WG(IGAUSS)*BK2YX*SFN*JACOB 
BROWY(M,2*IN)=BROWY(M,2*IN)+WG(IGAUSS)*BK2YY*SFN*JACOB 
END DO Each Gauss point in turn 

I 

! TERMS INVOLVING LOGARITHMIC QUADRATURE. 

Each Gauss point in turn: DO IGAUSS=1,NGAUSS 

IT=2 

IC=4 

CALL JACOBI(M,IGAUSS,IT,IC) 

SFN=SFL(IC,IN,IGAUSS) 
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DZDE=2. 

BROWX(M,2*IN-1)=BROWX(M,2*IN-1)+WGL(IGAUSS)*CL*SFN*JACOB*DZDE 
BROWY(M,2*IN)=BROWY(M,2*IN)+WGL(IGAUSS)*CL*SFN*JACOB*DZDE 
END DO Each Gauss point in turn 
END IF 
END IF 

END DO Each element node in turn 

RETURN 

END SUBROUTINE INTKER 
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P not at the current node of the element containing Q 

If P and Q are not in the same element, or P is not at the current node of the element containing Q 
(that is, I=£J) then normal Gaussian quadrature (Appendix A) over the whole element is carried out to 
evaluate the integrals. The required shape functions at the Gauss points have already been computed in 
subprogram SHAPE and stored in array SF. The Jacobians at the Gauss points are found in subprogram 
JACOBI, and the kernel functions in subprogram KERNEL. The last of these subprograms stores the 
kernel functions in variables AKXX, AKXY, AKYX, AKYY and BKXX, BKXY, BKYX, BKYY, respectively. 
These correspond to T xx , T xy , T yx , T yy and U xx , U xy , U yx , U yy in the equations. For each Gauss point 
the products of Gaussian weighting factor (WG), kernel functions, shape function (SFN) and Jacobian 
(JACOB) are then added to either arrays AROWX, AROWY, BROWX or BROWY to complete the 
numerical integration process. 

P at the current node of the element containing Q 

If P and Q are in the same element and P is at the current node of the element (that is, I=J) then, 
because of the singular nature of the second kernel function, some logarithmic Gaussian quadrature is 
required, as described in Section 5.5.2. The first kernel functions do not have to be evaluated because the 
relevant coefficients at the diagonal of the [A] matrix have already been found indirectly by summing 
the off-diagonal terms. 

If P is the first node of the element the integration process follows Equations 5.82 to 5.96. The integrals 
of the second kernel function are split into singular and non-singular parts. To the latter are applied 
Gaussian quadrature in the usual way, the non-singular parts of the second kernel functions being defined 
in subprogram KERN2 and returned in variables BK2XX, BK2XY, BK2YX and BK2YY. The singular 
parts are integrated using logarithmic Gaussian quadrature (Appendix A) applied to the whole element, 
with the origin of the modified intrinsic co-ordinate r\ being at the first node. The derivative of with 
respect to r/ , |d<f/d 77 |, is stored in variable DZDE. 

If P is the second node of the element the integration process follows Equations 5.97 to 5.109. The 
integrals of the second kernel functions are again split into singular and non-singular parts. To the latter 
are applied Gaussian quadrature in the usual way, the non-singular parts of the second kernel function 
being defined in subprogram KERN2 and returned in variables BK2XX, BK2XY, BK2YX and BK2YY. 
The singular parts are integrated using logarithmic Gaussian quadrature applied to the two halves of the 
element separately, in each case with the origin of the modified intrinsic co-ordinate 77 at the second node. 

If P is the third node of the element the integration process follows Equations 5.110 to 5.121, and is 
very similar to the case of P at the first node of the element. 

6.1.9 Subprogram for computing the Jacobian at a point on an element 

Subprogram JACOBI is identical to the subprogram described in Section 4.1.10, with the exception that 
the USE statement is for module SHAREDDATA 2 EQ rather than SHAREDDATA 2 PQ. 
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SUBROUTINE JACOBI(M,IGAUSS,IT,IC) 

i 

! SUBPROGRAM TO EVALUATE THE JACOBIAN AND THE COMPONENTS OF THE UNIT 
! NORMAL VECTOR AT A PARTICULAR GAUSS POINT. 

! M INDICATES THE ELEMENT NUMBER. 

! IGAUSS INDICATES THE GAUSS POINT NUMBER. 

! IT INDICATES THE TYPE OF QUADRATURE. 

! IT—1 - NORMAL GAUSSIAN QUADRATURE. 

! IT=2 - LOGARITHMIC GAUSSIAN QUADRATURE. 

! IC INDICATES THE CASE NUMBER FOR LOGARITHMIC GAUSSIAN QUADRATURE, 

! AS DEFINED IN SUBROUTINE SHAPE. 

i 

USE SHAREDDATA2EQ 

i 

! CALCULATE THE DERIVATIVES OF THE GLOBAL COORDINATES WITH RESPECT TO 
! THE LOCAL INTRINSIC COORDINATE. 

DX=0 . 

DY=0 . 

Each element node in turn: DO IN=1,3 
IF(IT == 1) SFND=SD(IN,IGAUSS) 

IF(IT == 2) SFND=SDL(IC,IN,IGAUSS) 

J=NODE(M,IN) 

DX=DX+SFND*XNODE(J) 

DY=DY+SFND*YNODE(J) 

END DO Each element node in turn 

i 

! COMPONENTS OF THE LOCAL OUTWARD NORMAL VECTOR AT THE GAUSS POINT. 
UNGX=DY 
UNGY=-DX 

I 

! JACOBIAN OF THE COORDINATE TRANSFORMATION. 

JACOB=SQRT(UNGX**2+UNGY**2) 

i 

! SCALE THE VECTOR COMPONENTS TO GIVE THE UNIT NORMAL VECTOR. 

UNGX=UNGX/JACOB 
UNGY=UNGY/JACOB 

I 

RETURN 

END SUBROUTINE JACOBI 
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By products of the calculation for a Jacobian (JACOB) are the components of the outward normal vector 
to the boundary at the chosen Gauss point. These are converted to components of the unit normal and 
stored in variables UNGX and UNGY. 

6.1.10 Subprogram for computing the kernel functions when P is not at the current node of the 
element containing Q 

Subprogram KERNEL computes the kernel functions at a particular Gauss point Q. The co-ordinates 
XQ and YQ of the point are first found using Equations 5.69 and 5.70, followed by the components r x 
and r y (RX and RY) of the radius vector of length r (R) joining P and Q. The rate of change of radius 
r with the normal n to the boundary (DRDN) is found as the cosine of the angle between them, or the 
scalar product of the two unit vectors. The first kernel functions are computed using Equations 5.33, 
5.36, 5.38, 5.40, and the second kernels from Equations 5.26 to 5.29. 


SUBROUTINE KERNEL(XP,YP,M,IGAUSS,UNX,UNY) 


SUBPROGRAM TO EVALUATE THE KERNELS WHEN P IS NOT THE CURRENT 
ELEMENT NODE. 

XP, YP INDICATE THE GLOBAL COORDINATES OF POINT P. 

M INDICATES THE ELEMENT NUMBER. 

IGAUSS INDICATES THE NUMBER OF THE GAUSS POINT, Q. 
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j 

USE SHAREDDATA2EQ 

| 

! COORDINATES OF GAUSS POINT Q. 

XQ=0 . 

YQ=0 . 

Each element node in turn: DO IN=1,3 
SFN=SF(IN,IGAUSS) 

J=NODE(M,IN) 

XQ=XQ+SFN*XNODE(J) 

YQ=YQ+SFN*YNODE(J) 

END DO Each element node in turn 

j 

! COMPONENTS AND MAGNITUDE OF THE RADIUS VECTOR FROM P TO Q. 

RX=XQ-XP 
RY=YQ-YP 
RSQ=RX**2+RY**2 
R=SQRT(RSQ) 

RX=RX/R 

RY=RY/R 

| 

! RATE OF CHANGE OF R WITH THE NORMAL TO THE BOUNDARY AT Q. 
DRDN=RX*UNX+RY*UNY 

j 

! PARAMETERS IN THE KERNEL FUNCTIONS. 

Cl=-1./(4.*PI*R) 

C2=l.-NU 
C3=2.*(l.+NU) 

C4=(1.+NU)**2/(4.*PI*E) 

C5=(3.-NU)*AL0G(1./R)/(l.+NU) 

j 

! EVALUATE THE KERNELS. 

AKXX=C1*(C2+C3*RX*RX)*DRDN 
TERM1=C3*RX*RY*DRDN 
TERM2=C2*(RX*UNY-RY*UNX) 

AKXY=C1*(TERM1-TERM2) 

AKYX=C1*(TERM1+TERM2) 

AKYY=C1*(C2+C3*RY*RY)*DRDN 
BKXX=C4 *(C5+RX*RX) 

BKXY=C4 *RX*RY 
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BKYX=BKXY 
BKYY=C4 *(C5+RY*RY) 

i 

RETURN 

END SUBROUTINE KERNEL 

6.1.11 Subprogram for computing the second kernel function when P is at the current node of 
the element containing Q 

Subprogram KERN2 computes, for a particular Gauss point Q , the non-singular components of the 
second kernel functions when force point is at the current node of the element containing Q. In other 
words, all but the ln^, ln^— anc ^ l n (“) terms in Equations 5.92 to 5.96, 5.108 and 5.109, 5.120 and 
5.121. These are stored in BK2XX, BK2XY, BK2YX and BK2YY in the subprogram. The value of function 
R (<f) (RFN in the subprogram) is defined by either Equation 5.91, 5.106 or 5.119, according to whether 
P is at the first, second or third node of the element. The relevant value of the intrinsic co-ordinate 
(variable ZETA in the program) is obtained from array ZG for normal Gaussian quadrature. 

SUBROUTINE KERN2(M,IN,IGAUSS) 

i 

! SUBPROGRAM TO EVALUATE THE NON-SINGULAR LOGARITHMIC TERM IN THE 
! SECOND KERNEL WHEN P IS THE CURRENT ELEMENT NODE. 

! M INDICATES THE ELEMENT NUMBER. 

! IN INDICATES THE NUMBER OF THE ELEMENT NODE FORMING P. 

! IGAUSS INDICATES THE GAUSS POINT NUMBER. 

i 

USE SHAREDDATA2EQ 

i 

! COORDINATES OF GAUSS POINT Q. 

XQ=0 . 

YQ=0 . 

Each element node in turn: DO 10—1,3 
SFN=SF(IC,IGAUSS) 

J=NODE(M,IC) 

XQ=XQ+SFN*XNODE(J) 

YQ=YQ+SFN*YNODE(J) 

END DO Each element node in turn 

I 

! ELEMENT NODE NUMBERS. 

Il=NODE(M,1) 

I2=NODE(M,2) 

I3=NODE(M,3) 

i 
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! EVALUATE THE INTRINSIC COORDINATE. 

ZETA=ZG(IGAUSS) 

| 

! P AT THE FIRST NODE OF THE ELEMENT. 

IF(IN == 1) THEN 

XCOMP=(ZETA-2.)*XNODE(II)+2.*(1.-ZETA)*XNODE(12)+ZETA*XNODE(13) 
YCOMP=(ZETA-2.)*YNODE(II)+2.*(1.-ZETA)*YNODE(12)+ZETA*YNODE(13) 
RFN=SQRT(XCOMP**2+YCOMP**2) 

1 = 11 
END IF 

| 

! P AT THE SECOND NODE OF THE ELEMENT. 

IF(IN == 2) THEN 

XCOMP=-0.5*(ZETA-1.)*XNODE(II)+ZETA*XNODE(12) 

& -0.5*(ZETA+1.)*XNODE(13) 

YCOMP=-0.5*(ZETA-1.)*YNODE(II)+ZETA*YNODE(12) 

& -0.5*(ZETA+1.)*YNODE(13) 

RFN=SQRT(XCOMP**2+YCOMP**2) 

1 = 12 
END IF 
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j 

! P AT THE THIRD NODE OF THE ELEMENT. 

IF(IN == 3) THEN 

XCOMP=-ZETA*XNODE(II)+2.*(1.+ZETA)*XNODE(12)-(ZETA+2.)*XNODE(13) 
YCOMP=-ZETA*YNODE(II)+2.*(1.+ZETA)*YNODE(12)-(ZETA+2.)*YNODE(13) 
RFN=SQRT(XCOMP**2+YCOMP**2) 

1 = 13 
END IF 

| 

! COMPONENTS AND MAGNITUDE OF THE RADIUS VECTOR FROM P TO Q. 

XP=XNODE(I) 

YP=YNODE(I) 

RX=XQ-XP 
RY=YQ-YP 
RSQ=RX**2+RY**2 
R=SQRT(RSQ) 

RX=RX/R 

RY=RY/R 

j 

! PARAMETERS IN THE KERNEL FUNCTIONS. 

C4=(1.+NU)**2/(4.*PI*E) 

C5=(3.-NU)*ALOG(1./RFN)/(1.+NU) 

| 

! EVALUATE THE KERNELS. 

BK2XX=C4 *(C5+RX*RX) 

BK2XY=C4 *RX*RY 
BK2YX=BK2XY 
BK2YY=C4*(C5+RY*RY) 

| 

RETURN 

END SUBROUTINE KERN2 

6.1.12 Results output subprogram 

Subprogram OUTPUT organises and writes out the primary results of the computation. Firstly, 
displacements and tractions are stored in their proper arrays. 

SUBROUTINE OUTPUT 

| 

! SUBPROGRAM TO WRITE OUT THE NODAL POINT VALUES OF DISPLACEMENTS 
! AND ELEMENT STRESSES AND COMPUTE THE FORCES ON THE BOUNDARY SEGMENTS. 

I 
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USE SHAREDDATA2EQ 

j 

! ARRANGE FOR U, V AND TX, TY TO CONTAIN THE NODAL DISPLACEMENTS 
! AND TRACTIONS, RESPECTIVELY. 

Each node in turn: DO 1=1,NNP 
IF(IBCN(I) == 1) THEN 

TX(I)=UV(2*1-1) 

TY(I)=UV(2*I) 

END IF 

IF(IBCN(I) == 2) THEN 
U(I)=UV(2*1-1) 

V(I)=UV(2*1) 

END IF 

END DO Each node in turn 

I 

! HEADING FOR OUTPUT OF NODAL DISPLACEMENTS AND TRACTIONS. 

WRITE(6,61) 

61 FORMAT(/ "NODAL POINT DISPLACEMENTS AND TRACTIONS" // 

& " I U V TX TY") 

E=ESTORE 

Each node in turn: DO 1=1,NNP 

| 

! REMOVE THE SCALING. 

U(I)=U(I)*MAXL 
V(I)=V(I)*MAXL 
TX(I)=TX(I)*ESTORE 
TY(I)=TY(I)*ESTORE 

| 

! OUTPUT THE NODAL VALUES OF DISPLACEMENTS AND TRACTIONS. 

WRITE(6,62) I,U(I) ,V(I) ,TX(I) ,TY (I) 

62 FORMAT(14,4E15.6) 

END DO Each node in turn 

I 

! HEADING FOR OUTPUT OF ELEMENT STRESSES. 

WRITE(6,63) 

63 FORMAT(/ "STRESSES AT THE NODES OF THE ELEMENTS" // 

& " M IN SIGNN SIGSS SIGSN", 

& " SIGE") 


97 


Download free eBooks at bookboon.com 


Boundary Element Methods for Engineers: 
Part II: Plane Elastic Problems 


Quadratic BoundaryElementProgramforPlaneElastic Problems 


! NORMAL AND SHEAR STRESSES AT THE NODES OF THE ELEMENTS. 

Each element in turn: DO M=1,NEL 
Each element node in turn: DO IN=1,3 
IF(IBCE(M) /= 2) THEN 

J=NODE(M,IN) 

TMX (M, IN) =TX (J) 

TMY (M, IN) =TY (J) 

SIGNN (M, IN) =TMX (M, IN) *UNMX (M, IN) +TMY (M, IN) *UNMY (M, IN) 

SIGSN (M, IN) =-TMX (M, IN) *UNMY (M, IN) +TMY (M, IN) *UNMX (M, IN) 

END IF 

| 

! DIRECT STRESS ALONG THE BOUNDARY SURFACE. 

IC1=N0DE(M,1) 

IC2=N0DE(M,2) 

IC3=N0DE(M,3) 

IGAUSS=IN+10 

DUDZ=SD(1,IGAUSS)*U(IC1)+SD(2,IGAUSS)*U(IC2)+SD(3,IGAUSS)*U(IC3) 
DVDZ=SD(1,IGAUSS)*V(IC1)+SD(2,IGAUSS)*V(IC2)+SD(3,IGAUSS)*V(IC3) 
IT=1 
IC=1 
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CALL JACOBI(M,10+IN,IT,IC) 

ESS=(-DUDZ*UNMY(M,IN)+DVDZ*UNMX(M,IN))/(JACOB*MAXL) 

SIGSS(M,IN)=E*ESS+NU*SIGNN(M,IN) 

i 

! VON MISES EQUIVALENT STRESS. 

SIGE(M,IN)=SQRT(SIGNN(M,IN)**2+SIGSS(M,IN)**2 
& -SIGNN(M,IN)*SIGSS(M,IN)+3.*SIGSN(M,IN)**2) 

i 

! OUTPUT THE STRESSES AT THE NODES OF THE ELEMENTS. 

WRITE(6,64) M,IN,SIGNN(M,IN),SIGSS(M,IN),SIGSN(M,IN), 

& SIGE(M,IN) 

64 FORMAT(214,4E15.6) 

END DO Each element node in turn 
END DO Each element in turn 

i 

! COMPUTE THE FORCES ON THE BOUNDARY SEGMENTS. 

Each segment in turn: DO ISEG=1,NSEGTOT 
FXSEG(ISEG)=0. 

FYSEG(ISEG)=0. 

END DO Each segment in turn 
Each element in turn: DO M=1,NEL 
ISEG=ISEGELEM(M) 

i 

! APPLY GAUSSIAN QUADRATURE. 

FXELEM=0. 

FYELEM=0. 

Each element node in turn: DO IN=1,3 

Each Gauss point in turn: DO IGAUSS=1,NGAUSS 

SFN=SF(IN,IGAUSS) 

IT-1 

IC=1 

CALL JACOBI(M,IGAUSS,IT,IC) 

FXELEM=FXELEM+WG(IGAUSS)*SFN*JACOB*TMX(M,IN)*MAXL 
FYELEM=FYELEM+WG(IGAUSS)*SFN*JACOB*TMY(M,IN)*MAXL 
END DO Each Gauss point in turn 
END DO Each element node in turn 

I 

! ACCUMULATE THE FORCES ON THE SEGMENT. 

FXSEG(ISEG)=FXSEG(ISEG)+FXELEM 
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FYSEG(ISEG)=FYSEG(ISEG)+FYELEM 
END DO Each element in turn 

i 

! OUTPUT THE SEGMENT FORCE RESULTS. 

WRITE(6,65) (ISEG,FXSEG(ISEG),FYSEG(ISEG),ISEG=1,NSEGTOT) 

65 FORMAT(/ "FORCES ACTING ON THE BOUNDARY SEGMENTS" 

& // "SEGMENT FX FY" / (I5,2E14.5)) 

i 

RETURN 

END SUBROUTINE OUTPUT 

If the boundary condition at a particular point is of the first type (prescribed displacements) the solution 
array UV contains the computed tractions, which are stored in arrays TX and TY. If the boundary 
condition is of the second type then displacements were the unknowns and are now stored in arrays 
U and V. Then all displacements and tractions have the scaling removed, and nodal point values of 
displacements and tractions are written out. 

Stresses at each of the three nodes of each element are then computed, unless they were given as input 
data. In many cases this will provide duplicate information for end nodes of adjacent elements, but not 
if such nodes are at corners of the domain where tractions are discontinuous. Normal and shear stresses 
are first found using Equations 5.60 and 5.61, followed by the direct stresses along the boundary from 
Equations 5.62 and 5.63. The von Mises equivalent stress is found from Equation 5.64, and all four 
stresses are written out. 

In many problems of practical engineering interest, it is useful to have the global co-ordinate components 
of the forces applied to each of the boundary segments. These can found by integrating the tractions 
element by element (forces stored in FXELEM and FYELEM), summing to give total segment forces, 
FXSEG and FYSEG, which are then written out. 

6.1.13 Subprogram for solving the linear equations 

Both the Gaussian elimination algorithm and subprogram ELIMIN for solving the linear equations are 
described in Appendix B. The matrix [A*] extended to include the right hand side vector [b ](in Equation 
5.57) is entered in array A, and the solutions returned in array X (UV in the calling main program). 

6.2 Some Test Problems for BEM2EQ 

Before attempting to solve real problems for which the solutions are not known, it is essential to test 
any computer program as extensively as possible on problems for which exact analytical solutions are 
available for comparison. Test problems serve not only to verify that the program is working correctly, 
but also to examine the level of discretisation (the fineness of the boundary element mesh in this case) 
necessary to obtain accurate results. 
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6.2.1 Simple tension 

One of the most basic test problems available is simple tension. Figure 6.1 shows a thin strip of material 
which has dimensions of 4 units by 2 units. Both the left and right hand edges of the strip (AD and 
BC) are subject to uniform normal tensile stresses of magnitude 3 units. The top and bottom edges (DC 
and AB) are free of applied stresses. Youngs modulus for the material of the strip is 2000 units, and its 
Poissons ratio is 0.3. Under plane stress conditions, with a uniaxial stress of a xx = 3 the direct strain in 
the x direction is e xx = ^ = 0.0015 and in they direction is e yy = —ve xx = —0.00045. Hence the 
length of the strip increases by 0.0015 X 4 = 0.006 and the width reduces by 0.00045 X 2 = 0.0009. 
The total forces acting on the left hand and right hand edges are both 6 (stress times edge length, the 
strip thickness being treated as unity). 
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Although the strip is in force equilibrium, the problem cannot yet be solved using boundary elements 
because the boundary conditions are given purely in terms of stresses, with no displacements prescribed. 
This difficulty can be overcome by applying suitable point constraints. The constraints adopted here are 
to fix the midpoint E of side AB in both the v and y directions, while constraining midpoint F of side 
DC not to move in the v direction (but to remain free to move in the y direction). The point constraints 
prevent rigid body movement of the strip as a whole, without giving rise to any point forces at either 
point E or point F. 

To obtain a quadratic boundary element solution to the problem, the boundary has first to be divided 
into segments. Four segments along the four edges of the domain are appropriate, because they are 
straight lines between corners and the boundary condition over each one is uniform. The ends of the 
four segments are the corner points labelled A, B, C and D (in that order) in Figure 6.1. The choice of 
starting point is arbitrary, but then the direction of numbering must keep the domain to the left. The 
data file DATA to solve the problem using one element (three nodes) per side of the domain is as follows 


SIMPLE TENSION 


STRESS 

(Plane stress selected) 

2000. 0.3 

(Youngs modulus and Poissons ratio) 

1 

(Number of boundaries) 

4 

(Number of segments on the boundary) 

0. 0. 4. 0. 4. 2. 0. 2. 

(Co-ordinates of the segment ends) 

0 . 1 1 . 

(Radii of curvature, 

0 . 11 . 

number of elements 

0 . 11 . 

and element length ratio 

0 . 11 . 

for each of the 4 segments) 

0 2 3 

(Two prescribed stresses, three point constraints) 

2 3. 0. 

(Normal stress = 3 on segment 2) 

4 3. 0. 

(Normal stress = 3 on segment 4) 

2 1 

(Node 2 constrained in v direction) 

2 2 

(Node 2 constrained in y direction) 

6 1 

(Node 6 constrained in v direction) 


The annotations on the right have been added for explanation, and are not part of the data file. Note that 
the prescribed zero values of stresses (on DC and AB) do not have to be entered explicitly. 


102 


Download free eBooks at bookboon.com 


Boundary Element Methods for Engineers: 
Part II: Plane Elastic Problems 


Quadratic BoundaryElementProgramforPlaneElastic Problems 


The mesh data output file MESHRES is as follows 


GEOMETRIC DATA FOR THE MESH 
NUMBER OF ELEMENTS = 4 
NUMBER OF NODAL POINTS = 8 
COORDINATES OF THE NODAL POINTS 


I 

X 

Y 

I 

X 

Y 

1 

0.00 0 0E + 0 0 

0.0000E+00 

2 

0.20 0 0E + 01 

0.0000E+00 

3 

0.4000E + 01 

0.0000E + 00 

4 

0.40 0 0E + 01 

0.1000E+01 

5 

0.4000E + 01 

0.2000E + 01 

6 

0.2000E+01 

0.2000E+01 

7 

0.0000E + 00 

0.2000E + 01 

8 

0.00 0 0E + 0 0 

0.1000E+01 


ELEMENT 

NODE 

NUMBERS 





M 

ND1 

ND2 

ND3 

M 

ND1 

ND2 

ND3 

1 

1 

2 

3 

2 

3 

4 

5 

3 

5 

6 

7 

4 

7 

8 

1 


With one element per segment, 4 elements and 8 nodes are created. Both elements and nodes are 
numbered as for the segments, anticlockwise from the origin, as shown in Figure 6.2. Element numbers 
are shown circled. Element end nodes are shown as small solid circles, mid-side nodes as open circles. 



The RESULTS data file output by the program is 

QUADRATIC BOUNDARY ELEMENT SOLUTION FOR TWO DIMENSIONAL ELASTIC PROBLEM 
SIMPLE TENSION 

UNDER PLANE STRESS CONDITIONS 

YOUNG'S MODULUS = 0.2000E+04 POISSON'S RATIO = 0.300 
PRESCRIBED STRESS BOUNDARY CONDITIONS 

SIGNN = 0.3000E+01 SIGSN = 0.0000E+00 ON ELEMENTS 2 TO 2 
SIGNN = 0.3000E+01 SIGSN = 0.0000E+00 ON ELEMENTS 4 TO 4 
NODE 2 CONSTRAINED WITH U = 0 
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NODE 2 CONSTRAINED WITH V = 0 
NODE 6 CONSTRAINED WITH U = 0 
NODAL POINT DISPLACEMENTS AND TRACTIONS 


I U 

1 -0.2 99965E-02 

2 0.000000E+00 

3 0.2 99965E-02 

4 0.300182E-02 

5 0.2 99965E-02 

6 0.000000E+00 

7 -0.299965E-02 

8 -0.300182E-02 


V 

-0.101365E-05 
0.000000E+00 
-0.101577E-05 
-0.450233E-03 
-0.899451E-03 
-0.900459E-03 
-0.899450E-03 
-0.450231E-03 


TX 

-0.300000E+01 
0.000000E+00 
0.300000E+01 
0.300000E+01 
0.150000E+01 
0.000000E+00 
-0.300000E+01 
-0.300000E+01 


TY 

0.000000E+00 
0.000000E+00 
-0.799680E-06 
0.000000E+00 
-0.399840E-06 
0.000000E+00 
0.000000E+00 
0.000000E+00 


STRESSES AT THE NODES OF THE ELEMENTS 


M IN SIGNN 

1 1 0.000000E+00 
1 2 0.000000E+00 

1 3 0.000000E+00 

2 1 0.300000E+01 
2 2 0.300000E+01 

2 3 0.300000E+01 

3 1 0.000000E+00 
3 2 0.000000E+00 

3 3 0.000000E+00 

4 1 0.300000E+01 
4 2 0.300000E+01 
4 3 0.300000E+01 


SIGSS 

0.299965E+01 
0.2 99965E+01 
0.2 99965E+01 
0.156770E-02 
0.15 6531E-02 
0.15 62 87E-02 
0.2 99965E+01 
0.299965E+01 
0.2 99965E+01 
0.15 622 9E-02 
0.15 6415E-02 
0.15 65 95E-02 


SIGSN 

0.000000E+00 
0.000000E+00 
0.000000E+00 
0.000000E+00 
0.000000E+00 
0.000000E+00 
0.000000E+00 
0.000000E+00 
0.000000E+00 
0.000000E+00 
0.000000E+00 
0.000000E+00 


SIGE 

0.299965E+01 
0.299965E+01 
0.299965E+01 
0.299922E+01 
0.299922E+01 
0.299922E+01 
0.299965E+01 
0.299965E+01 
0.299965E+01 
0.299922E+01 
0.299922E+01 
0.299922E+01 


FORCES ACTING ON THE BOUNDARY SEGMENTS 


SEGMENT 

FX 

FY 

1 

0.00000E+00 

0.00000E+00 

2 

0.60000E+01 

-0.53312E-06 

3 

0.00000E+00 

0.00000E+00 

4 

-0.60000E+01 

0.00000E+00 
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The prescribed stress conditions are a normal value of 3 and shear value of zero over element 2 (nodes 3, 
4 and 5) on the right end edge (segment 2), and a normal value of 3 and shear value of zero over element 
4 (nodes 7, 8 and 1) on the left hand edge (segment 4). Segments (and elements) 1 and 3 defining the 
top and bottom edges had no boundary condition prescribed by the DATA file, so zero stress conditions 
are assumed. 


The computed displacements show that at nodes 3, 4 and 5 (on the right hand edge) the displacements 
in the direction are 0.0029996, 0.0030018 and 0.0029996, respectively, while at nodes 7, 8 and 1 (on the 
left hand edge) they are -0.0029996, -0.0030018 and -0.0029996. The extensions of the strip between 
node pairs 3 and 1, 4 and 8 and 5 and 7 are therefore 0.0059992, 0.0060036 and 0.0059992, respectively. 
These compare very favourably with the exact value of 0.006. Note the slight differences between element 
end nodes and element midside nodes, which were also experienced in the quadratic element analysis of 
potential problems. The displacement in the y direction of node 6 (which is vertically above constrained 
node 2) is -0.00090046, which compares with the exact value of -0.0009. 


The value of o X x (SIGSS) computed at all the nodal points on the top and bottom edges is 2.9996, 
compared with the exact value of 3. Also, small tensile values of o yy (SIGSS) of about 0.00156 are 
computed at all the nodes on the right and left hand side edges, where they should be zero. The forces on 
the left and right hand edges are exact at 6.0000, which is to be expected because they were prescribed 
in the boundary conditions. 
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The results are highly accurate for this very simple problem. Using just one quadratic element and three 
nodes on each of the four edges of the boundary gives results which are almost correct to four significant 
figures. For this very simple problem it is not necessary to refine the mesh because the accuracy of the 
results is already adequate for most practical engineering purposes. 

Using displacement boundary conditions 

An alternative approach to the simple tension problem shown in Figure 6.1 is to apply displacement 
boundary conditions to one end of the strip, say edge AD. Unfortunately, this normally changes the 
nature of the problem because contraction along AD is prevented. The only exception is if the value of 
Poissons ratio is set to zero, rendering the problem one-dimensional. 

The RESULTS data file output for this case is 


QUADRATIC BOUNDARY ELEMENT SOLUTION FOR TWO DIMENSIONAL ELASTIC PROBLEM 
SIMPLE TENSION 

UNDER PLANE STRESS CONDITIONS 

YOUNG'S MODULUS = 0.2000E+04 POISSON'S RATIO = 0.000 

PRESCRIBED DISPLACEMENT BOUNDARY CONDITIONS 
U = 0.0000E + 00 V = 0.0000E + 00 ON ELEMENTS 4 TO 4 
PRESCRIBED STRESS BOUNDARY CONDITIONS 

SIGNN = 0.3000E + 01 SIGSN = 0.0000E + 00 ON ELEMENTS 2 TO 2 

NODAL POINT DISPLACEMENTS AND TRACTIONS 


I U V 

1 0.000000E+00 0.000000E+00 

2 0.300082E-02 0.376321E-07 

3 0.60007 6E-02 -0.847336E-06 

4 0.6002 60E-02 -0.344766E-08 

5 0.60007 6E-02 0.840583E-06 

6 0.300083E-02 -0.407261E-07 

7 0.000000E+00 

8 0.000000E+00 


TX 

-0.300536E+01 
0.000000E+00 


TY 

0.254296E-03 
0.000000E+00 


0.300000E+01 -0.799680E-06 
0.300000E+01 0.000000E+00 

0.150000E+01 -0.3 9984 0E-0 6 
0.000000E+00 0.000000E+00 

0.000000E+00 -0.300536E+01 -0.251318E-03 
0.000000E+00 -0.299728E+01 -0.181283E-06 


STRESSES AT THE NODES OF THE ELEMENTS 


M IN SIGNN SIGSS 

1 1 0.000000E+00 0.300127E+01 
1 2 0.000000E+00 0.300038E+01 

1 3 0.000000E+00 0.299949E+01 

2 1 0.300000E+01 0.1687 64E-02 
2 2 0.300000E+01 0.168792E-02 


SIGSN SIGE 

0.000000E+00 0.300127E+01 
0.000000E+00 0.300038E+01 
0.000000E+00 0.299949E+01 
0.000000E+00 0.299916E+01 
0.000000E+00 0.299916E+01 
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2 3 0.300000E+01 

3 1 0.000000E+00 
3 2 0.000000E+00 

3 3 0.000000E+00 

4 1 0.300536E+01 
4 2 0.2 99728E+01 
4 3 0.300536E+01 


0.168820E-02 
0.299949E+01 
0.300038E+01 
0.300127E+01 
0.000000E+00 
0.000000E+00 
0.000000E+00 


0.000000E+00 
0.000000E+00 
0.000000E+00 
0.000000E+00 
0.251318E-03 
0.181283E-06 
-0.254296E-03 


0.299916E+01 
0.2 9994 9E+01 
0.300038E+01 
0.300127E+01 
0.300536E+01 
0.2 99728E+01 
0.300536E+01 


FORCES ACTING ON THE BOUNDARY SEGMENTS 


SEGMENT 

FX 

FY 

1 

0.00000E+00 

0.00000E+00 

2 

0.60000E+01 

-0.53312E-06 

3 

0.00000E+00 

0.00000E+00 

4 

-0.60000E+01 

0.7507 9E-0 6 


The extensions of the strip computed at nodes 3, 4 and 5 on edge BC are 0.0060007, 0.0060026 and 
0.0060007 (exact 0.006), the o xx stresses at nodes on the top and bottom edges of the strip are 3.0013, 
3.0004 and 2.9995 (exact 3), and the force on edge AD is exact at 6.0000. The accuracy remains good 
when displacements are prescribed rather than stresses. 


6.2.2 Pure shear 

In order to explore shear loading, Figure 6.3 shows the cross-section of a rectangular block of material 
which is long in the direction normal to the cross-section. The material concerned again has a Youngs 
modulus of 2000 and a Poissons ratio of 0.3. The base of the block, AB, is rigidly constrained with 
prescribed zero values of displacements. The other edges have shear stresses of magnitude 3 applied in 
such a way as to give rise to a state of pure shear. In other words, on sides BC and DA a sn = +3, and 
on side CD a sn = — 3, bearing in mind that co-ordinate S is measured along the boundary keeping the 
domain always to the left. 


C 

A 

x — 3 

B 

Figure 6.3 Block under pure shear 


D 


r = 3 
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The shear modulus for the material is 


E _ 2000 

2(1 + v) _ 2(1 + 0.3) 


769.23 


and the shear strain due to a shear stress of 3 units is 3/769.23 = 0.0039. With a block height of 2, the 
displacement of top surface DC is u = 2 X 0.0039 = 0.0078. 


The RESULTS data file output for this case is 


QUADRATIC BOUNDARY ELEMENT SOLUTION FOR TWO DIMENSIONAL ELASTIC PROBLEM 
PURE SHEAR 

UNDER PLANE STRAIN CONDITIONS 

YOUNG'S MODULUS = 0.2000E+04 POISSON'S RATIO = 0.300 
PRESCRIBED DISPLACEMENT BOUNDARY CONDITIONS 
U = 0.0000E+00 V = 0.0000E+00 ON ELEMENTS 1 TO 1 
PRESCRIBED STRESS BOUNDARY CONDITIONS 


SIGNN = 

0.0000E+00 

SIGSN = 

0.3000E+01 

ON 

ELEMENTS 

2 

TO 

2 

SIGNN = 

0.0000E+00 

SIGSN = 

-0.3000E+01 

ON 

ELEMENTS 

3 

TO 

3 

SIGNN = 

0.0000E+00 

SIGSN = 

0.3000E+01 

ON 

ELEMENTS 

4 

TO 

4 
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NODAL POINT DISPLACEMENTS AND TRACTIONS 


I 

U 

V 

TX 

TY 

1 

0.000000E+00 

0.000000E+00 

-0.300535E+01 

-0.184857E-02 

2 

0.000000E+00 

0.000000E+00 

-0.299735E+01 

0.206969E-06 

3 

0.000000E+00 

0.000000E+00 

-0.300535E+01 

0.184803E-02 

4 

0.390155E-02 

0.144327E-05 

0.000000E+00 

0.300000E+01 

5 

0.7 800 91E-02 

0.117580E-06 

0.150000E+01 

0.150000E+01 

6 

0.780417E-02 

0.283357E-09 

0.300000E+01 

0.000000E+00 

7 

0.780091E-02 

-0.118 911E-0 6 

0.150000E+01 

-0.150000E+01 

8 

0.390155E-02 

-0.144314E-05 

0.000000E+00 

-0.300000E+01 


STRESSES 

AT 

THE NODES OF 

THE ELEMENTS 



M 

IN 

SIGNN 

SIGSS 

SIGSN 

SIGE 

1 

1 

0.184857E-02 

0.792242E-03 

-0.300535E+01 

0.520542E+01 

1 

2 - 

0.206969E-06 

-0.887008E-07 

-0.299735E+01 

0.519156E+01 

1 

3 - 

0.184803E-02 

-0.792012E-03 

-0.300535E+01 

0.520542E+01 

2 

1 

0.000000E+00 

0.621710E-02 

0.300000E+01 

0.519616E+01 

2 

2 

0.000000E+00 

0.12 92 09E-03 

0.300000E+01 

0.519615E+01 

2 

3 

0.000000E+00 

-0.595411E-02 

0.300000E+01 

0.519616E+01 

3 

1 

0.000000E+00 

-0.716504E-02 

-0.300000E+01 

0.519616E+01 

3 

2 

0.000000E+00 

0.000000E+00 

-0.300000E+01 

0.519615E+01 

3 

3 

0.000000E+00 

0.716504E-02 

-0.300000E+01 

0.519616E+01 

4 

1 

0.000000E+00 

0.595146E-02 

0.300000E+01 

0.519616E+01 

4 

2 

0.000000E+00 

-0.130672E-03 

0.300000E+01 

0.519615E+01 

4 

3 

0.000000E+00 

-0.621280E-02 

0.300000E+01 

0.519616E+01 


FORCES ACTING ON THE BOUNDARY SEGMENTS 


SEGMENT 

FX 

FY 

1 

-0.12000E+02 

0.19397E-06 

2 

0.53312E-06 

0.60000E+01 

3 

0.12000E+02 

-0.26656E-06 

4 

0.00000E+00 

-0.60000E+01 


The computed displacements at nodes 5, 6 and 7 on the top surface DC are 0.0078009, 0.0078042 and 
0.0078009 (exact 0.0078). The computed shear stresses at nodes 1, 2 and 3 on the bottom surface AB 
are -3.00535, -2.9973 and -3.0053 (exact 3), and the total shear force on this surface is exact at -12.000. 
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It is worth noting that had the shear stresses on sides BC and DA been omitted in Figure 6.3, although the 
block might appear to be simply sheared, the actual state of stress would have been much more complex. 
In particular, at corners A and B there is a conflict between finite shear stress along the constrained base 
and zero shear stress on the side of the block. Equation 1.1 requires shear stresses to be complementary. 

6.2.3 Steel component 

Figure 6.4 shows a thin steel component with elastic properties of Youngs modulus 207 GN/m 2 , Poissons 
ratio 0.30. The bottom edge is subject to a uniform normal stress of 100 MN/m 2 , and the top edges 
to uniform normal stresses of 75 MN/m 2 (giving force equilibrium in the vertical direction). The top 
edges are also subject to outwardly directed shear stresses of 30 MN/m 2 , giving force equilibrium in the 
horizontal direction. 



This problem involves somewhat more complex geometry and varying stresses, and is therefore likely to 
require rather more mesh refinement than the very simple problems treated so far. The object is to find 
the direct stress in the surface at point E (mid way between points D and F on the semi-circular arc), 
and the maximum von Mises equivalent stress anywhere on the boundary (and hence anywhere in the 
component). There is no exact analytical solution for comparison. 

The boundary of the domain is conveniently divided into six segments: AB, BC, CD, DEF, FG and GA. 
For convenience, the same number of boundary elements are used on each segment. Displacement 
constraints are required to prevent rigid body movement: for example, point A fixed in both directions, 
and point B constrained not the move in the direction. Table 6.1 shows the effects of varying the number 
of nodes per segment from 4 to 40. Less than 4 gives poor geometric representation of the semi-circle 
DEF (with each element then representing an arc of more than 45°). 
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Elements 

a ss at E 

max.Og 

per segment 

MN/m 2 

MN/m 2 

4 

-111.3 

152.2 

(140.2) 

6 

-119.9 

151.5 

(148.7) 

10 

-123.8 

153.7 

(153.7) 

20 

-125.1 

152.9 

(152.9) 

40 

-125.3 

152.7 

(152.7) 


Table 6.1 Stress results for varying levels of mesh refinement 


In each case point E, which is on the vertical line of symmetry for the problem, is located at the common 
node of two elements, and the value of & S s obtained from each element is the same (to at least four 
significant figures). On the other hand, the maximum value of equivalent stress, which occurs in the 
regions of both points H and I in Figure 6.4, does not necessarily occur at a node at all. This means that 
a reasonable degree of mesh refinement is needed to capture a good approximation to the maximum. In 
each case the maximum is detected at an element end node, and the two values of equivalent stress from 
the two elements meeting at that node are not necessarily the same. In the table the values in brackets 
are those from the adjacent element. 

There are two criteria to keep in mind when testing for convergence with mesh refinement. Firstly, 
whether a value of stress changes with increasing number of elements. Secondly, whether stress values 
at an element end node (in this case the maximum equivalent stress) are the same for both elements 
meeting at that node. With these in mind, the second criterion in this problem is satisfied at 10 elements 
per segment or more, whereas the first is only satisfied (to three significant figures) at 20 elements per 
segment or more. The computed direct stress in the surface at point E is 125 MN/m 2 (compressive) and 
the maximum von Mises equivalent stress 153 MN/m 2 (actually the tensile direct stresses in the surfaces 
in the regions of points H and I). 

6.2.4 A thick-walled cylinder test problem 

The next test is designed to demonstrate the ability of the program to solve problems with curved 
boundaries and those with more than one boundary. Figure 6.5 shows the cross-section of a long thick- 
walled circular cylinder with internal and external radii of iq — 4 and r 2 = 8. The cylinder is subject to 
an internal pressure of p = 3 and no pressure on the outside. The Youngs modulus and Poissons ratio 
of the material of the cylinder are 2000 and 0.3, respectively. Cylinders are widely used in engineering 
practice to contain fluids under pressure: from pipes for relatively low pressures to thick-walled vessels 
for high-pressure chemical reactors. 
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J 



The analytical solution to the thick-walled cylinder problem, under either plane stress or plane strain 
conditions, gives a general result for the radial and hoop stresses as 

B B 

°rr = A - Ogg = A + ^ (6.5) 
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where r is radial distance from the centre of the cylinder, and and are constants. Stress o rr is the direct 
radial stress in the r direction while hoop stress Oqq is, as its name suggests, the direct stress in the 
circumferential direction at any given radial position (for example, along the inner or outer surfaces of 
the cylinder). The values of the constants depend on the boundary conditions applied to the cylinder. 
In the present problem, these are 

o rr — —p at r = r 1 and o rr — 0 at r - r 2 

Hence A = f ; and B = where K = — 

(K 2 - 1) (X 2 -l) r*i 

and the hoop stresses at the inner and outer cylinder surfaces are 

p{K 2 + 1) , 2p , 

°oo = iK 2 _ x) at r = and o ee = at r = r 2 ( 6 . 8 ) 

In this problem K — 2 and p = 3, and these hoop stresses are 5 and 2 , respectively. 

These results for stresses do not depend on the values of either Youngs modulus or Poissons ratio. Radial 
displacement, u r , can be found from the hoop strain 


eee = ( a h ~ W r ) (6.9) 

At the inner cylinder surface u r = (5 + 3v*) (6.10) 

E* 

and at the outer surface u r = —7 ( 6 . 11 ) 

t 

From Equations 5.7 

E* = —Tt = 2197.80 and v* = -X- = 0.428571 (6.12) 

(l-v 2 ) (1—v) 


and the displacements at the inner and outer cylinder surfaces are 

u r = 0.011440 and u r = 0.007280 


( 6 . 6 ) 

(6.7) 
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Bearing in mind that program BEM2EQ cannot accept boundary segment specifications with angles 
greater than , the simplest mesh to represent the geometry shown in Figure 6.5 requires four segments, 
two on the outer boundary with end points arbitrarily chosen at (8, 0) and (-8, 0), and two on the inner 
boundary with end points at (4, 0) and (-4, 0). These end points are the points F, G, H and I in Figure 
6.5. The origin for the global Cartesian co-ordinates is the centre of the cylinder. The following is an 
annotated DATA input file for this problem, for the case of four elements (9 nodes) per semi-circular 
segment, two per quadrant. This means that each element represents an arc of by means of a parabolic 
(quadratic) curve passing through both the ends and midpoint of the arc. 

THICK-WALLED CYLINDER 


STRAIN 

(Plane strain selected) 

2000. 0.3 

(Youngs modulus and Poissons ratio) 

2 

(Number of boundaries) 

2 

(Number of segments on first boundary) 

o 

00 

1 

o 

00 

(Co-ordinates of the segment ends) 

8. 4 1. 

(Radii of curvature, numbers of elements 

00 

I — 1 

and length ratios) 

2 

(Number of segments on second boundary) 

o 

i 

o 

(Co-ordinates of the segment ends) 

\—i 

i 

(Radii of curvature, numbers of elements 

-4. 4 1. 

and length ratios) 

0 2 3 

(2 prescribed stresses and 3 point constraints) 

3 -3. 0. 

(Normal stress = -3 on segment 3) 

4 -3. 0. 

(Normal stress = -3 on segment 4) 

1 1 

(Node 1, point F, constrained in x direction) 

1 2 

(Node 1, point F, constrained in y direction) 

9 2 

(Node 9, point G, constrained in y direction) 


As in the case of the simple tension problem of Section 6.2.1, sufficient displacements constraints must 
be prescribed to prevent rigid body movement of the whole cylinder. Those chosen are no movement 
in either the or direction for point F, and no movement in the direction for point G. In terms of node 
numbers, for the present mesh of 4 elements per segment these are 1 and 9. 
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Table 6.2 shows the results obtained for meshes ranging from one element per quadrant up to four 
elements per quadrant. For the computed stresses, two values are shown, the first for element end nodes, 
and the second for midside nodes. The computed displacements are for the highest points on surfaces, 
J and K in the figure, in the direction. 


Elements 

At inner surface 

At outer surface 

per quadrant 

G ee 

u r 

G ee 

u r 

Exact 

5 

0.011440 

2 

0.007280 

1 

4.8985 

0.011029 

1.9653 

0.007030 


4.7743 


1.9313 


2 

4.9940 

0.011417 

1.9967 

0.007213 


4.9875 


1.9965 


3 

5.0006 

0.011437 

1.9996 

0.007279 


4.9984 


1.9996 


4 

5.0014 

0.011440 

2.0001 

0.007280 


5.0001 


2.0001 



Table 6.2 Results for the thick-walled cylinder problem 
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Since the problem is axi-symmetric it is to be expected that the computed stresses at all the nodes on the 
inner boundary would be identical, and similarly for all the nodes on the outer boundary. This proves to 
be the case, but for expected fluctuations between the values at element end nodes and midside nodes. 
At four elements per quadrant the results are almost correct to four significant figures, which is adequate 
for most practical purposes. In this case each element covers a circular arc of 22.5°, which demonstrates 
the ability of quadratic elements to accurately model curved boundaries. Thirty two elements having 64 
nodes are sufficient to model this multiply-connected solution domain. 

6.3 An Example: Confined Compression of a Rubber Block 

A practical problem for solution by program BEM2EQ is the compression of a rubber block which is 
confined between rigid surfaces. It has important engineering applications in the bearings for bridges. 
Bridges are often supported on rubber pads. These are relatively flexible in shear, allowing the bridge to 
expand or contact in length with changes in ambient temperature and the load it is carrying. But they 
are much stiffer in compression. The reason for this is that rubber, which is much more flexible than 
metals, is incompressible (or very nearly incompressible), having a Poissons ratio of 0.5. It is worth 
noting that many finite element methods are incapable of analysing incompressible material problems 
under plane strain conditions. 

•i 


Figure 6.6 

Figure 6.6 shows a rubber block sandwiched between rigid surfaces, to which it is bonded. The height 
and width of the block are H and W, respectively. Consider the case where W/H = 10.. The dimension 
of the block in the direction normal to the plane shown is large, so that a state of plane strain may be 
assumed. With the bottom surface AB fixed in both global co-ordinate directions, the top surface is to 
be lowered by 0.1% of the thickness H. The (maximum) compressive stress at point E is to be computed, 
together with the total compressive force on the surface AB, both to be expressed as multiples of the 
stress and force which would exist if the block were in a state of simple compression. 

If the Youngs modulus is taken as 1000, a 0.1% simple compressive strain would give a stress of 
0.001 X 1000 = 1. If block width is taken as 1, then the total forces on the upper and lower surfaces 
would also be 1 under the same conditions. So in the boundary element model it is convenient to 
take W = 1, H = 0.1 and an upper surface vertical displacement of v — —0.0001 (0.1% of 0.1). The 
computed stresses and forces will be the required multiples. 
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Rubber block between rigid surfaces 
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The boundary of the domain is conveniently divided into four segments: AB, BC, CD and DA. With 
a wide disparity in the lengths of the segments, it is no longer appropriate to use the same number of 
elements on each segment. Best results are normally obtained when the lengths of adjacent elements are 
reasonably similar (not differing by more than a factor of about 2 or 3). In this problem it is appropriate 
to have ten times as many elements on AB and CD as on BC and DA. This means that the coarsest mesh 
which can be used has 10 elements each on AB and CD, and only one each on BC and DA. This is also 
desirable from the point of view that the domain is slender and it is important that opposite sides are 
not too close together in relation to element sizes (for the reasons discussed in Section 4.2.1 for potential 
problems). With 10 elements on each of the horizontal edges, their distance apart is equal to the common 
lengths of the elements. Mesh refinement (keeping the 10:1 element number ratio between the horizontal 
and vertical segments) will only improve the situation. 

Table 6.3 shows the effects of varying the number of elements on the horizontal segments from 10 to 50 
(with the numbers of elements on the vertical segments always one tenth of these). 


Elements 

per segment 

Oyy at E 

Total force 

on AB 

10 

-50.43 

-34.51 

20 

-50.15 

-34.03 

30 

-50.18 

-33.95 

40 

-50.19 

-33.92 

50 

-50.19 

-33.91 


Table 6.3 Stress and force results for varying levels of mesh refinement 


There is no change in the maximum stress of 50.2 beyond 20 elements per horizontal segment and 
the total force has settled to 33.9 by 30 elements per segment. As multiples of the values expected for 
simple compression these are large numbers, and would have been much higher if larger ratios had 
been considered. The explanation is that, because the rubber material is assumed to be completely 
incompressible, the reduction in the height of the block caused by bringing the confining surfaces closer 
together can only be accommodated by material being forced out horizontally at the side edges, BC 
and DA. Because there are no displacements at the corners A, B, C and D the rubber is forced out into 
curved bulges. In this case the maximum horizontal displacements at the centres of BC and DA are each 
computed as 0.00072, 7.2 times the relative displacement of the confining surfaces. The driving forces 
required to produce this deformation are provided by the stress gradients in the horizontal direction 
moving away from the centre of the block. From a maximum of 50.2 at point E, both o xx and a (the 
state of stress in the block is hydrostatic, with these stresses equal) follow a roughly parabolic horizontal 
distribution to finish at close to one at sides BC and DA. The average stress is 33.9, numerically equal 
to the total applied force. 
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6.4 An Example: Stress Concentration at a Hole in a Flat Plate 

As a final example, consider the classic problem of the stress concentration created at a small circular 
hole at the centre of a thin flat square plate under tension. This is often used as a test of any numerical 
method of stress analysis. Figure 6.7 shows the geometry and loading. The radius of the hole is , and the 
width and height of the plate are 2: the ratio of hole diameter to plate width is . For present purposes this 
ratio is taken as 1/20. A uniform tensile stress is applied to both the top and bottom edges of the plate. 



Figure 6.7 Square plate with a small central hole 

There is no exact analytical solution for the case of a plate of finite width and height, although one does 
exist for an infinitely large plate. This gives a maximum (tensile) hoop stress at both points E and G of 
3a, together with a (compressive) hoop stress at both points F and H of -a. It also indicates that there 
are rapid changes in stresses and deformations close to the hole, which make this a relatively challenging 
application for any numerical method. The present problem is for a plate of finite width, although the 
width to diameter ratio is relatively high. The maximum stress concentration factor (ratio of local stress 
to remote applied stress) can be expected to be a little higher than 3. 

The outer boundary of the domain is conveniently divided into four segments: AB, BC, CD and DA, and 
the inner boundary into two: EHG and GFE. On the outer boundary the stresses and displacements do 
not vary much and only a few boundary elements are required: four per segment is sufficient. Rather 
more elements are likely to be required on the inner boundary. Point constraints are required to prevent 
rigid body movement: point A is constrained in both directions, while point B is constrained not to 
move in the y direction. 
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The DATA file for this problem with four elements per semi-circular segment (two per quadrant) is 


STRESS CONCENTRATION 

STRESS 

2000. 0.3 

2 

4 

-0.5 -0.5 0.5 -0.5 0.5 0.5 -0.5 0.5 
0. 4 1. 

0. 4 1. 

0. 4 1. 

0. 4 1. 

2 

0.025 0. -0.025 0. 

-0.025 4 1. 

-0.025 4 1. 

0 2 3 

1 1 . 0 . 

3 1. 0. 

1 1 
1 2 
9 2 

SIMPLY CLEVER 


(Plane stress selected) 

(Youngs modulus and Poissons ratio) 

(Number of boundaries) 

(Number of segments on first boundary) 
(Co-ordinates of the segment ends) 

(Radii of curvature, numbers of elements 
and length ratios) 

(Number of segments on second boundary) 
(Co-ordinates of the segment ends) 

(Radii of curvature, numbers of elements 
and length ratios) 

(2 prescribed stresses and 3 point constraints) 
(Normal stress = 1 on segment 1) 

(Normal stress = 1 on segment 3) 

(Node 1 constrained in x direction) 

(Node 1 constrained in y direction) 

(Node 9 constrained in y direction) 

SKODA 
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The overall dimensions of the plate are taken as 1, so that a=0.025 and b= 0.5. The Youngs modulus and 
Poissons ratio are arbitrarily chosen as 2000 and 0.3, respectively, although they do not affect the results 
for stress concentration. 


An edited version of the RESULTS data file output produced is 


QUADRATIC BOUNDARY ELEMENT SOLUTION FOR TWO DIMENSIONAL ELASTIC PROBLEM 

STRESS CONCENTRATION 

UNDER PLANE STRESS CONDITIONS 

YOUNG'S MODULUS = 0.2000E+04 POISSON'S RATIO = 0.300 

PRESCRIBED STRESS BOUNDARY CONDITIONS 

SIGNN = 0.1000E+01 SIGSN = 0.0000E+00 ON ELEMENTS 1 TO 4 

SIGNN = 0.1000E+01 SIGSN = 0.0000E+00 ON ELEMENTS 9 TO 12 

NODE 1 CONSTRAINED WITH U = 0 

NODE 1 CONSTRAINED WITH V = 0 

NODE 9 CONSTRAINED WITH V = 0 

STRESSES AT THE NODES OF THE ELEMENTS 


M 

IN 

SIGNN 

SIGSS 

SIGSN 

SIGE 

17 

1 

0.000000E+00 

0.301684E+01 

0.000000E+00 

0.301684E+01 

17 

2 

0.000000E+00 

0.242750E+01 

0.000000E+00 

0.242750E+01 

17 

3 

0.000000E+00 

0.105800E+01 

0.000000E+00 

0.105800E+01 

18 

1 

0.000000E+00 

0.945216E+00 

0.000000E+00 

0.945216E+00 

18 

2 

0.000000E+00 

-0.424362E+00 

0.000000E+00 

0.424362E+00 

18 

3 

0.000000E+00 

-0.10137 0E+01 

0.000000E+00 

0.101370E+01 

19 

1 

0.000000E+00 

-0.10137 0E+01 

0.000000E+00 

0.10137 0E+01 

19 

2 

0.000000E+00 

-0.424375E+00 

0.000000E+00 

0.424375E+00 

19 

3 

0.000000E+00 

0.945128E+00 

0.000000E+00 

0.945128E+00 

20 

1 

0.000000E+00 

0.105804E+01 

0.000000E+00 

0.105804E+01 

20 

2 

0.000000E+00 

0.242751E+01 

0.000000E+00 

0.242751E + 01 

20 

3 

0.000000E+00 

0.301677E+01 

0.000000E+00 

0.301677E+01 

21 

1 

0.000000E+00 

0.301691E+01 

0.000000E+00 

0.301691E+01 

21 

2 

0.000000E+00 

0.242751E+01 

0.000000E+00 

0.242751E + 01 

21 

3 

0.000000E+00 

0.105795E+01 

0.000000E+00 

0.105795E+01 

22 

1 

0.000000E+00 

0.945229E+00 

0.000000E+00 

0.94522 9E+00 

22 

2 

0.000000E+00 

-0.424373E+00 

0.000000E+00 

0.424373E+00 

22 

3 

0.000000E+00 

-0.101371E+01 

0.000000E+00 

0.101371E+01 

23 

1 

0.000000E+00 

-0.101371E+01 

0.000000E+00 

0.101371E+01 

23 

2 

0.000000E+00 

-0.424386E+00 

0.000000E+00 

0.424386E+00 

23 

3 

0.000000E+00 

0.945244E+00 

0.000000E+00 

0.945244E+00 

24 

1 

0.000000E+00 

0.105804E+01 

0.000000E+00 

0.105804E+01 

24 

2 

0.000000E+00 

0.242751E + 01 

0.000000E+00 

0.242751E+01 

24 

3 

0.000000E+00 

0.301683E+01 

0.000000E+00 

0.301683E+01 

FORCES 

l ACTING ON THE 

BOUNDARY SEGMENTS 
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SEGMENT 

FX 

FY 

1 

-0.56196E-07 

-0.10000E+01 

2 

0.00000E+00 

0.00000E+00 

3 

0.56196E-07 

0.10000E+01 

4 

0.00000E+00 

0.00000E+00 

5 

0.00000E+00 

0.00000E+00 

6 

0.00000E+00 

0.00000E+00 


Results for nodal point displacements and tractions have been omitted, together with stresses on the 
outer boundary (elements 1 to 16). The hoop stress (SIGSS) at point E is either that for the first node of 
element 17 (3.01684) or the third node of element 24 (3.01683). At point G it is either that for the third 
node of element 20 (3.01677) or the first node of element 21 (3.01691). The four values are consistent, 
and give an average to five significant figures of 3.0168. The hoop stress at point F is either that for the 
third node of element 18 (-1.01370) or the first node of element 19 (-1.01370). At point H it is either 
that for the third node of element 22 (-1.01371) or the first node of element 23 (-1.01371). The four 
values are consistent, and give an average to five significant figures of -1.01371. With only two elements 
per quadrant of the hole in the plate the stress values are very close to what is expected. 
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Elements 

per quadrant 

Gyy 

at E and G 

°xx 

at F and H 

2 

3.0168 

-1.0137 

3 

3.0206 

-1.0163 

4 

3.0211 

-1.0167 

6 

3.0215 

-1.0168 


Table 6.4 Stress concentration results for varying levels of mesh refinement at the hole 

Table 6.4 shows the effects of refining the mesh around the hole. Both stress values increase in magnitude 
a little and converge to 3.021 at the sides of the hole and -1.017 at the top and bottom of the hole. 
Converged results are achieved with only a modest number of boundary elements: at 6 elements per 
quadrant the total used is 40. With rapid variations of stresses and deformations near the hole, this is the 
type of problem where the ability to vary the element distributions along segments (to give small elements 
at the stress concentrations) might have been expected to offer significant benefits. The quality of the 
results from using elements of constant size suggests, however, that this is not necessary. The situation 
would have been different if symmetry had been invoked to model only one quarter of the plate (say, the 
first quadrant in the plane). Element size variation on the lines of symmetry would have been required 
to accommodate both small elements near the hole and much larger elements at the outer boundary. 

6.5 Discussion 

It is clear that quadratic boundary elements applied to elastic stress analysis problems offer good 
accuracy with relatively small numbers of elements. The scope of the program described in this chapter 
is deliberately limited, to keep the coding reasonably straightforward to understand. Many other features 
can be incorporated: for example, additional types of boundary conditions, displacement and stress 
calculations at internal points, and inclusion of body forces such as gravity and centrifugal effects. 

Another useful facility is the ability is to be able to divide a domain into several distinct regions, with 
elements along the boundaries between them. The different regions can be of different materials, having 
different elastic properties. This makes possible the treatment of not only multi-material problems, but 
also multiple bodies in contact with each other. Another class of problems that can be catered for is that 
of crack problems in linear elastic fracture mechanics. One of the boundary element end nodes can be 
located at the crack tip and the midside nodes in the adjacent elements shifted to force the required stress 
singularity at the crack tip (just as is done in finite element methods). If both faces of the crack have to 
be modelled, then the multi-region facility is very useful. This is because coincident nodal points in a 
single boundary element region are not permitted, and nearly coincident nodes give rise to inaccurate 
analysis. Coincident and nearly coincident nodes in different regions, on the two sides of the crack, 
present no difficulties. 
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Problems 

6.1 In the stress concentration problem illustrated in Figure 6.7 (with a/b = 1/2 0)> instead of the 
tensile loading applied to the top and bottom edges of the plate, the hole itself is subjected to 
a uniform internal pressure. Use program BEM2EQ to find the maximum stress concentration 
factor (relative to the applied pressure). Compare this value with the available analytical solutions. 

6.2 Consider the long thick-walled cylinder whose cross-section is shown in Figure 6.5, but with its 
outer surface rigidly contained. An internal pressure of p = 1 is applied. The Youngs modulus 
and Poissons ratio of the material are 1000 and 0.3, respectively. Use program BEM2EQ to find 
both the stresses in the cylinder at the inner and outer surfaces, and the radial displacement of 
the inner surface, and compare with analytical results. 

6.3 Figure 6.8 shows a thin rectangular plate with a small semi-circular notch at the centre of one 
horizontal edge. The radius of the notch is l/20 th of the vertical height of the plate. The plate 
is subjected to uniform tensile stresses on its vertical edges. Use program BEM2EQ to find the 
maximum stress concentration at the notch. Take Poissons ratio to be 0.4. 



Figure 6.8 Notched plate in tension 
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6.4 Figure 6.9 shows a thin square plate with a square hole at its centre. The corners of the hole are 
rounded, with radii of curvature 10% of the hole dimensions, 5% of the overall plate dimensions. 
Use program BEM2EQ to find the maximum stress concentration at the hole. Take Poissons 
ratio to be 0.25. 



6.5 A cast iron pedestal is trapezoidal in cross-section, with a horizontal base 0.3 m wide, horizontal 
top 0.1 m wide and height 0.1 m, and is symmetrical about the vertical centre line. Plane strain 
conditions may be assumed. The Youngs modulus is 70 GN/m 2 and Poissons ratio 0.25. The top 
surface of the pedestal is uniformly loaded in compression and the base rests on a rigid surface. 
The purpose of the pedestal is to spread the load over the supporting surface. Use program 
BEM2EQ to find the ratio between the maximum compressive stress applied to the supporting 
surface and the stress applied to the top of the pedestal. The surface is rough enough to prevent 
any slipping. Is a coefficient of friction of 0.5 large enough to achieve this? 
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6.6 Figure 6.10 shows a steel tension test piece, which is 3 mm thick. The central section is 20 
mm wide, increasing to 30 mm at the ends where it is gripped. The transition in width is via 
quadrant fillets with radii of 5 mm. Youngs modulus and Poissons ratio for steel may be taken 
as 207 GN/m 2 and 0.3, respectively. The ends of the test piece are subjected to uniform stresses 
of 60 MN/m 2 . Use program BEM2EQ to find the stiffness of the test piece and the maximum 
stress induced. Find also the maximum and minimum axial stresses over the central 30 mm, 
the region where an extensometer is to be attached and where it is intended that the axial stress 
is uniform. 


a 

30 

1 


20 ► 5 R 


20 , 




20 

J. _/ 


r 


1 -v 


100 


60 MN/m 2 


All dimensions in mm 


Figure 6.10 Tension test piece 
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6.7 If the edges of the plate in Figure 6.7 (with alb - 1/20) are subject to pure shear rather than 
simple tension, use program BEM2EQ to find the maximum stress concentration at the hole, 
expressed as a multiple of the applied shear stress. 

6.8 Consider the stress concentration problem shown in Figure 6.7 not with a circular hole at the 
centre but with an elliptical hole. The major axis of the ellipse lies along the x axis and is l/20 th 
of the plate width, while the minor lies along the axis and is half the major axis. Use program 
BEM2EQ to find the maximum stress concentration factor. Take Poissons ratio to be 0.3. 

6.9 Consider the stress concentration problem shown in Figure 6.7 not with a single circular hole at 
the centre but with two holes of equal size (alb = 1/20), side by side. Their centres are on the x 
axis and equidistant from the centre of the plate. The distance between the centres is four times 
their common radius. Use program BEM2EQ to find the stress concentration factors at the sides 
of the holes. Take Poissons ratio to be 0.3. Consider both (a) when the applied uniform tensile 
stresses are in the y direction as shown, at right angles to the line of hole centres, and (b) when 
they are in the x direction, parallel to the line of centres. 

6.10 Figure 6.11 shows the cross-section of a long circular steel cylinder, internal and external radii 
100 and 250 mm, respectively. Running the full length of the cylinder is a straight groove which 
is semi-circular in shape with a radius of 5 mm. Youngs modulus and Poissons ratio for steel 
may be taken as 207 GN/m 2 and 0.3, respectively. Use program BEM2EQ to find the maximum 
von Mises equivalent stress in the cylinder when an internal pressure of 120 MN/m 2 is applied. 


All dimensions in mm 
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6.11 The long thick-walled internally pressurised cylinder described in Section 6.2.4 has been 
incorrectly manufactured so that it is slightly elliptical in cross-section, with major axis 1% 
greater than minor axis. Use program BEM2EQ to find the ratio between the maximum hoop 
stress (at the internal surface) and the applied pressure, and compare this to the value for the 
intended circular cylinder. 

6.12 Part of the suspension for a motor vehicle involves a torsional spring formed by a circular 
cylindrical rubber bush which is 75 mm long and has inner and outer diameters of 20 mm and 
50 mm, respectively. The outer surface is bonded to a rigid housing, and the inner surface to a 
metal shaft. The rubber has a Youngs modulus of 80 MN/m 2 and a Poissons ratio of 0.5. Use 
program BEM2EQ to find the torque which must be applied to the shaft to rotate it by 0.1 radians. 
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7 Further Applications 

In Chapter 1 many continuum mechanics problems in a wide range of engineering subjects were 
categorised according to their governing differential equations into either the harmonic or biharmonic 
types. Problems governed by harmonic equations are, in engineering terms, more often referred to as 
potential problems, such as heat conduction, ideal fluid flow, porous medium flow, torsion and at least 
one form of fluid flow. Part I of the book was devoted to such potential problems, and boundary element 
methods were developed. Part II has been concerned with biharmonic problems in elastic stress analysis, 
particularly two-dimensional plane strain and plane stress. The purpose of this final chapter is to review 
briefly some further applications of boundary element methods. 

7.1 Axi-symmetric Problems 

There are many problems in engineering, both harmonic and biharmonic, which are symmetrical in terms 
of both geometry and boundary conditions about some axis. For example, the thick-walled cylinder test 
problems considered in Sections 3.2.2, 4.2.2 and 6.2.4 were symmetrical about the axis of the cylinder. In 
each case a cross-section of the cylinder at right angles to the axis was treated as the solution domain. It 
would also have been possible to treat as the solution domain a radial cross-sectional plane through the 
wall of the cylinder, which also passed through the cylinder axis. Boundary elements could be used to 
model the boundary of the domain. Interestingly, if a solid rather than hollow axi-symmetric component 
is modelled, so that the axis of symmetry forms part of the boundary of the solution domain, then it is 
not necessary to have boundary elements on the axis. The level of mathematics required to formulate 
axi-symmetric analyses is, however, rather more challenging than for two-dimensional problems. 

7.2 Higher-Order Boundary Elements 

It is of course possible to use boundary elements with shape functions of orders higher than quadratic. 
For example, some work has been done with cubic elements, including some which provide a smooth 
transition (continuity of slope) between one element and the next. The high quality of the results that 
can be obtained using quadratic elements, however, means that they are generally the elements of choice. 

7.3 Three-dimensional Problems 

Three-dimensional problems are amenable to solution by a boundary element approach. Again, only 
the boundary of the relevant solution domain is modelled. The surface is covered with two-dimensional 
boundary elements. The equivalent in three-dimensional analysis of the constant element for two- 
dimensional problems is the constant triangle or quadrilateral. Over either such element the potential 
and potential gradient normal to the surface (in a potential problem) or the displacements and tractions 
(in a stress analysis problem) are assumed to be constant over the element, and equal to the value at 
the node at the middle of the element. For potential problems there remains only one unknown at each 
node. For stress analysis there are now three, typically either displacements or tractions in each of the 
three global co-ordinate directions, and three equations are required at each node. 
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Progressing from constant to quadratic elements, the latter can again be either triangular or quadrilateral 
in shape, though not necessarily flat or with straight edges (in the same way that quadratic line elements 
are not necessarily straight lines). Isoparametric elements are generally preferred, with the geometric 
shape being defined in terms of quadratic functions of the three spatial co-ordinates. Triangular elements 
typically have six nodes, quadrilateral eight, one at each corner and one at the centre of each side. 
Quadrilateral elements are often used to model most of a surface, with triangles being useful to deal 
with local geometric features, particularly where a mesh refinement is required to suit a particular 
problem. Indeed, there are situations where a mesh refinement cannot be effected without using at least 
a few triangles. 

While three-dimensional meshes of (curved) triangles and quadrilaterals are significantly more difficult 
to create than two-dimensional meshes of line elements, they still only involve the discretisation of the 
boundary rather than the interior of a domain as well. Using boundary elements, the dimensionality 
of the computational problem is reduced by one. Nevertheless, the amount of computation involved in 
a three-dimensional boundary element analysis is very substantial, particularly because the coefficient 
matrix in the final set of algebraic equations is essentially fully populated with non-zero values. 
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7.4 Non-Linear Problems 

All the problems so far considered are linear in the sense that they involve the solution of sets of linear 
algebraic equations. Non-linearities can be introduced by either geometric, material property or boundary 
condition effects. Geometric non-linearities arise, for example, in problems involving solid media in 
which the strains are sufficiently large to significantly affect the shape of the solution domain. Such 
problems can be solved by an incremental approach in which the non-linear analysis is replaced by a 
series of linear analyses for progressively increasing external loads, after each of which the boundary 
element mesh geometry is recomputed. Non-linear behaviour due to boundary conditions arises in, for 
example, heat transfer problems with radiation boundary conditions, where the rate of heat transfer to 
or from a boundary depends on the fourth power of the absolute temperature there, rather than the first 
power in a convection boundary condition. In elastic stress analysis problems, contact between two or 
more bodies gives rise to non-linear behaviour: increased contact force between them changes the region 
of contact and hence modifies the boundary conditions. Again incremental approaches are possible. 

Examples of material non-linearities include non-Newtonian fluids and non-linearly elastic or elastoplastic 
solids, whose properties are functions of the local state of deformation or stress. For example, the viscosity 
of a fluid may be a function of the local strain rates (velocity gradients), or the elastic modulus of a solid 
may be a function of the local state of stress or strain. The form of dependence in each case would be 
determined from experimental data. It is also possible for material properties, such as viscosity, thermal 
conductivity, or elastic modulus to be significantly dependent on temperature, and therefore not known 
in advance. Such problems cannot be solved by boundary element method which discretise only the 
boundary: some internal discretisation and analysis is required adequately to treat the non-linearity. At 
least some of the benefits of the boundary element approach may thereby be lost. 

7.5 Comparison with Other Methods 

Other methods for solving engineering continuum mechanics problems include finite element, finite 
difference and finite volume methods. They are what are often referred to as domain methods: the entire 
two- or three-dimensional domain of the problem has to be discretised, rather than just the boundary. 
The principal advantage of boundary element methods is the effective reduction in dimensionality of the 
numerical problem. At least for linear situations, a two-dimensional physical problem requires only a 
one-dimensional mesh of boundary line elements, while a three-dimensional physical problem requires a 
two-dimensional mesh of typically triangular or quadrilateral elements. Benefits of this simplification can 
be seen in the problems addressed in this book: meshes of line elements for two-dimensional problems 
are easily generated using only a small number of geometric parameters, leading to very concise data 
input files. Analyses of even quite complicated problems can be set up quickly and easily. 
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Although it has not been demonstrated in this book, boundary element meshes, particularly of quadratic 
elements, are not merely domain method meshes with the internal points removed. For a given level of 
boundary discretisation, boundary element methods normally give substantially greater accuracy. But, 
although far fewer algebraic equations have to be solved, they have to be treated as being fully populated 
with non-zero coefficients. The computational benefits, although still present, are less clearcut. 

Features of boundary element methods which are perhaps deterrents to their wider use by engineers 
are the mathematics involved. Not only are they arguably more complex than for domain methods, but 
also tend to be of unfamiliar types: for example, integral equations, singular fundamental solutions, and 
analyses often expressed in tensor notation. An aim of this book has been to try to simplify and clarify the 
mathematics involved, even when this is sometimes at the expense of more bulky and less mathematically 
elegant equations. A further aim has been to demonstrate that boundary element methods are worthy 
of serious consideration for engineering analyses, and are particularly straightforward to use. 
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Appendix A: Gaussian Quadrature 

The aim of this appendix is to explain the Gaussian quadrature method of numerical integration, 
particularly as applied to boundary elements. 

When a function cannot be integrated analytically, some form of numerical integration or quadrature 
must be used. Two of the simplest methods are the familiar trapezium and Simpsons rules. Given a 
function /(x), values of the function are obtained for a number of values of the independent variable x 
over the required range of integration, and an approximate value of the integral obtained from them. An 
important feature (and limitation) of these simple methods is that the values of the independent variable 
at which the function is sampled are often taken at constant intervals, and certainly little attention is 
paid to the optimum choice of sampling points. 

If there is no restriction on the positions at which the function can be evaluated, then there is a class 
of numerical integration procedures referred to as Gaussian quadrature, which is to be preferred. The 
advantage over the trapezium and Simpsons rules is that greater accuracy for a given number of function 
evaluations, and therefore of overall computational effort, can be obtained. 



Figure A.1 Typical function to be integrated numerically 


Figure A.l shows a typical function to be integrated: the area between the curve /(<f) and the <f axis 
between the limits <f = — 1 and = +1 is required. It is convenient to transform the actual independent 
variable, such as x, into a dimensionless variable with these limits. This is the reason for introducing the 
local intrinsic co-ordinate for quadratic boundary elements. Assuming that the function is evaluated 
at n points, known as the Gauss points, the required integral can be expressed approximately as 

l\V(Od^Z2 =1 co p /(f p ) (A.1) 
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where a) p the are weighting factors. Simpsons rule (defined in Equation 3.49) is actually a particular 
form of this result in which, for a single application of the rule 

n = 3 

?l - - 1 ^2=0 (3 = +1 (A.2) 

1 4 

( 0 1 = 0 ) 3 = - "2=3 

Given the freedom to choose the <f p , however, the values of these and the corresponding weighting 
factors can be optimised for the integration of polynomial functions to give better accuracy. These are 
now well established and are available for a wide range of n> the number of Gauss points (for example, 
Stroud & Secrest 1966). 

For any numerical integration, choosing the number of Gauss points involves a compromise between 
accuracy and the amount of computation involved. In boundary element methods the number of points 
is normally at least n = 4. In this book a higher value of n = 8 is preferred, for which the values of 
and 0) are 


^ = 0.9602898564 

= -f 2 = 0.7966664774 

zz -f 3 zz 0.5255324099 

= -<f 4 = 0.1834346424 


o) 1 = o) 8 = 0.1012285362 
0) 2 = (o 7 = 0.2223810344 
a) 3 = a> 6 = 0.3137066458 
a) 4 = o) 5 = 0.3626837833 


(A.3) 


to ten decimal places. 

Gaussian quadrature formulae have also been worked out for functions which involve particular forms of 
functions as multipliers, including cases of multipliers which are singular within the range of integration. 
One such case which is very relevant to boundary element methods is 

Jo 1 In (i) f0i)dr] » £p=i 6%f(ri*) (A.4) 
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the numerical values for n = 8 being 


= 0.013320244 
r)2 = 0.079750429 
7 ) 3 = 0.197871029 
7)* 4 = 0.354153994 


o)| = 0.164416605 
c o* 2 = 0.237525610 
&>3 = 0.226841984 
u>% = 0.175754079 


7)1 = 0.529458575 
7)1 = 0.701814530 
7)* 7 = 0.849379320 
7)g = 0.953326450 


rug = 0.112924030 
col = 0.057872211 
co* 7 = 0.020979074 
col = 0.003686407 


(A.5) 


Reference 

Stroud, A.H. & Secrest, D. 1966, Gaussian Quadrature Formulas , Prentice-Hall. 
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Appendix B: Gaussian Elimination 

The aim of this appendix is to explain Gaussian elimination, which is a direct method for solving sets 
of simultaneous linear algebraic equations of the general form 

a 11 X 1 + a 12 X 2 + . + CLin x n — 

a 21 x l + a 22 x 2 + . + a 2n x n ~ b2 

. (B.l) 


CL n l x l T CL n2 X 2 T . T ^ nn x n b n 

where x 1 ,x 2 , m ■ ■ a ,x n are the unknowns, and the coefficients cty and b t are all known constants. This 
set of equations can be expressed in matrix form as follows 


' a ll 

a 21 

a 12 

a 22 

CLin 

a 2n 


-x t - 

*2 


[foil 

fo 2 

-^nl 

a n2 

CL nn - 


-X n - 


-bn- 


[A][x] = [b] 


(B.2) 


(B.3) 


The unknowns are successively eliminated by algebraic manipulation. The first equation can be used 
to eliminate x 1 from the remaining n — 1 equations. The modified second equation is then used to 
eliminate x 2 from the remaining n — 2 equations, and so on until the last equation contains only x n - 
Thus may be found, followed by all the other unknowns, by back substitution. Let the coefficients of [A] 
and [ b ] shown in Equations B.l and B.2 be given the notation a® and b^- After the nth elimination, 
the modified coefficients are 


a <r> = af-0<> 

- 0i,f) 

l IK 
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where 0 = CL ik /a kk -, i = k + l,k + 2, ...,n and j = k,k + 1, Note that the column vector [b] 
is treated just like a column of [/I ], and advantage can be taken of this fact to simplify the computer 
programming of the process. The final set of equations is 

a[ 1 1 ) x 1 + a^x 2 + .+ cS^x n = 

„(. 2) v + + „(2) v _ h (2) 

a 22 X 2 ' . ' a 2n X n ~ °2 

. (B.5) 


in) b (n) 

^nn u n 


Expressed in matrix terminology, the elimination process triangularises [A]. The unknowns are obtained 
in reverse order 


X =b (n) /a (n) 

A n U n / U nn 


(B.6) 



where i = n - l,n - 2,....,1 


A great many arithmetic operations are involved in solving large sets of equations by elimination. Any 
errors introduced, such as roundoff errors caused by representing numbers with only a finite number 
of significant figures, tend to be magnified and may become unacceptably large. Equations B.4 show 
that the elimination process involves many multiplications by the factors 0. In order to minimize the 
effects of roundoff errors, these factors should be made as small as possible, and certainly less than one. 
Therefore, the pivotal coefficient should be the largest coefficient in the leading column of the 
remaining submatrix 


i (k) 

l kk 


> 


(fc) 


a ik 


for i = k + 1 , k + 2 , ...,n 


(B.7) 


This condition also helps to avoid division by zero in Equations B.6, and can be achieved by a process 
known as partial pivoting. Immediately before each elimination, the leading column is searched for the 
largest coefficient. By interchanging equations, this can be made the pivotal coefficient to satisfy Equation 
B.7. The idea of partial pivoting can be extended to searching the whole of the remaining submatrix 
for the largest coefficient. Such complete pivoting involves interchanging both rows and columns, and 
is more difficult to program. Since it offers only modest advantages in terms of accuracy over partial 
pivoting, it is rarely used. 
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Another refinement which helps to minimize the effects of roundoff errors is to scale the equations to 
make their coefficients similar in magnitude. One way to do this is to normalize each equation so that 
the largest coefficient in each row of [A] is of magnitude one. Scaling can be particularly important 
when corresponding coefficients in different equations differ by several orders of magnitude. 

Any given set of equations may provide either duplicate (for example, two equations identical) or 
inconsistent information, and cannot be solved. In either case, the matrix of coefficients [A] is said to 
be singular, and its determinant is zero. Such a condition can be detected during partial pivoting if it 
is impossible to find a non-zero pivotal value at some stage of the elimination process or at the start of 
the back substitution process. What is more difficult to detect, however, is when a set of equations is 
nearly singular or ill-conditioned. This arises when, for example, two equations provide very nearly the 
same information about the unknowns. Alternatively, the effect of roundoff errors may be to make what 
is a singular set of equations apparently ill-conditioned, in that rather than a zero pivotal coefficient, a 
very small value is obtained. Indeed, the detection of a very small pivotal coefficient can be used as a 
convenient test for either a singular or very ill-conditioned set of equations, but without being able to 
distinguish between them. By very small pivotal coefficient is meant a value very small in relation to the 
magnitudes of the coefficients of the matrix at the start of the elimination process, which are of order 
unity if scaling has been applied. 

The following subprogram ELIMIN provides a Fortran implementation of Gaussian elimination. 

SUBROUTINE ELIMIN(A,X,NEQN,NROW,NCOL,IFLAG) 

I 

! SUBROUTINE FOR SOLVING SIMULTANEOUS LINEAR EQUATIONS BY GAUSSIAN 

! ELIMINATION WITH PARTIAL PIVOTING. 

i 

REAL :: A(NROW,NCOL),X(NROW) 

i 

! INITIALIZE ILL-CONDITIONING FLAG. 

IFLAG=0 

i 

! SCALE EACH EQUATION TO HAVE A MAXIMUM COEFFICIENT MAGNITUDE OF UNITY. 

JMAX=NEQN+1 

Each equation in turn: DO 1=1,NEQN 

AMAX=0. 

Search for maximum: DO J=1,NEQN 

ABSA=ABS(A(I,J)) 

IF(ABSA > AMAX) AMAX=ABSA 

END DO Search for maximum 
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Scale coefficients: DO J=1,JMAX 

A(I, J) =A(I, J) /AMAX 

END DO Scale coefficients 

i 

END DO Each equation in turn 

i 

! COMMENCE ELIMINATION PROCESS. 

Eliminate each variable in turn: DO K=1,NEQN-1 

i 

! SEARCH LEADING COLUMN OF THE COEFFICIENT MATRIX FROM THE DIAGONAL 
! DOWNWARDS FOR THE LARGEST VALUE. 

IMAX=K 

Search for largest value: DO I=K+1,NEQN 
IF(ABS(A(I,K)) > ABS(A(IMAX,K))) IMAX=I 
END DO Search for largest value 

i 

! IF NECESSARY, INTERCHANGE EQUATIONS TO MAKE THE LARGEST COEFFICIENT 
! BECOME THE PIVOTAL COEFFICIENT. 

IF(IMAX /= K) THEN 
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Interchange coefficients: DO J=K, JMAX 
ATEMP=A (K, J) 

A(K, J) =A (I MAX, J) 

A(IMAX,J)=ATEMP 
END DO Interchange coefficients 
END IF 

i 

! ELIMINATE X(K) FROM EQUATIONS (K+l) TO NEQN, FIRST TESTING FOR 

! EXCESSIVELY SMALL PIVOTAL COEFFICIENT (ASSOCIATED WITH A SINGULAR 

! OR VERY ILL-CONDITIONED MATRIX). 

IF(ABS(A(K, K)) < l.E-5) THEN 
IFLAG=1 
RETURN 
END IF 

Each of remaining equations: DO I=K+1,NEQN 
FACT=A(I,K)/A (K, K) 

Modify coefficients: DO J=K, JMAX 
A (I, J) =A(I, J)-FACT*A(K,J) 

END DO Modify coefficients 

END DO Each of remaining equations 

i 

END DO Eliminate each variable in turn 

i 

! SOLVE THE EQUATIONS BY BACK SUBSTITUTION, FIRST TESTING 

! FOR AN EXCESSIVELY SMALL LAST DIAGONAL COEFFICIENT. 

IF(ABS(A(NEQN,NEQN)) < l.E-5) THEN 
IFLAG=1 
RETURN 
END IF 

X(NEQN)=A(NEQN,JMAX)/A(NEQN,NEQN) 

Then each unknown in turn backwards: DO I=NEQN-1,1,-1 
SUM=A(I,JMAX) 

Sum products: DO J=I+1,NEQN 
SUM=SUM-A(I,J)*X(J) 

END DO Sum products 
X(I)=SUM/A(I,I) 

END DO Then each unknown in turn backwards 
RETURN 

END SUBROUTINE ELIMIN 
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The arguments of ELIMIN include the array of equation coefficients, A, which incorporates the vector 
of constants [b] as its last column, and the solution vector, X. The argument NEQN enters the number 
of equations to be solved, while IFLAG returns an integer flag value to the calling program to warn of 
a singular or very ill-conditioned coefficient matrix. 

The first action of the program is to scale the coefficients of the equations to give a maximum coefficient 
magnitude of unity in each row of [A ]. The elimination process is then started: each equation, with 
the exception of the last, is used in turn to eliminate the corresponding unknown from the subsequent 
equations. The current elimination number is given by K, and is equivalent to k in Equations B.4. Before 
performing the necessary eliminations with a particular equation, however, a search is made down the 
leading column of the remaining submatrix to find the coefficient of greatest magnitude, as defined 
by Equation B.7. The search technique locates the row number of the largest coefficient, IMAX, by 
first assuming that it corresponds to the coefficient on the diagonal, and only changing this if a larger 
coefficient is found. If the pivotal coefficient is not the largest, the relevant equations are interchanged 
by interchanging all their coefficients. In computing terms, this movement of data between storage 
locations is inefficient. While it could be avoided by keeping a record of the revised order in which the 
equations are to be considered, this makes the program significantly more difficult to understand. For 
present purposes, some computational efficiency is sacrificed in the interests of clarity. 

Despite the search for the largest coefficient, it is still possible for the pivotal coefficient to be extremely 
small or zero if the coefficient matrix is singular or very ill-conditioned. Bearing in mind that the equations 
were scaled initially, an appropriate test of magnitude is \a kk \ < 10 -5 . If this condition is satisfied at any 
stage of the elimination process, the problem is rejected. Rejection is indicated by setting the value of 
IFLAG to one, which may be detected by the calling program. If the partial pivoting is successful, however, 
the eliminations defined by Equations B.4 are performed, with the variable FACT being used to store 
the values of the factor 0. After testing the magnitude of the last diagonal coefficient (a^ n in Equations 
B.5), the back substitutions defined in Equations B.6 are performed to find the required solutions. 
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Appendix E: Matlab Version of 
Quadratic Boundary Element 
Program for Plane Elastic Problems 


The aim of this appendix is to present a Matlab version of the program BEM2EQ which is described in 
detail in its Fortran form in Chapter 6. Matlab is now very widely used by engineers for many purposes. 

The Matlab version is as far as possible a literal translation of the Fortran. This means that the detailed 
explanations provided in Chapter 6 are also applicable to the Matlab script. On the other hand, no 
attempt is made to take advantage of special built-in features of Matlab such as simplified handling of 
matrices and vectors, and the solution of sets of linear equations. 

Fortran is a compiler which converts the source code into machine code before computation starts. In 
contrast, Matlab is an interpreter which goes through the code line-by-line, translating and executing 
each line as it goes. Strictly speaking therefore, Matlab works with scripts rather than programs, although 
scripts are often loosely referred to as programs. A consequence of this difference between a compiler and 
an interpreter is that for large computationally intensive problems Fortran generally executes noticeably 
more quickly than Matlab. 
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The same function and variable names are used in the Fortran and Matlab programs: upper case names 
in Fortran become lower case in Matlab. The definitions of the names given in Some Program Variable 
Names at the beginnings of both Parts of the book are equally valid for both. DO loops in Fortran become 
“for” loops in Matlab, and functions are called in somewhat different ways. Also, STOP in Fortran is 
replaced by an error message to the computer screen in Matlab. The main difference, however, is to be 
seen in the details of the handling of input and output, although the same approach of inputting data 
from one pre-defined file and outputting to other files is adopted. In particular, the required input DATA 
files used in Fortran are unchanged, and the output files are effectively identical. 

The Matlab version of BEM2EQ, stored as the file BEM2EQ.m, is as follows. 

function BEM2EQ(varargin) 

o, 

o 

% PROGRAM FOR SOLVING TWO DIMENSIONAL ELASTIC STRESS ANALYSIS PROBLEMS 
% BY THE BOUNDARY ELEMENT METHOD USING QUADRATIC ELEMENTS. 

o 

o 

clear global; clear functions; 
global fid5 fid6 fid7 

o, 

o 

shareddata2eq 
fid5=f open ( ' DATA' , ' r ' ) ; 
fid6=f open ( ' RESULTS ' , ' w+ ' ) ; 
fid7=f open ( ' MESHRES ' , ' w+ ' ) ; 

o, 

o 

% DEFINE THE MAXIMUM PROBLEM SIZE PERMITTED BY THE ARRAY DIMENSIONS. 

maxnel=250; 

maxnnp=500; 

maxnb=10; 

maxnpc=2 0; 

o, 

o 

% INPUT THE PROBLEM TITLE AND TYPE, ALSO MATERIAL PROPERTIES, 
intitle; 

o 

o 

% INPUT AND GENERATE THE MESH DATA, 
meshq; 

o, 

o 

% OUTPUT THE MESH DATA, 
mshout; 

o 

o 
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% EVALUATE AND STORE VALUES OF THE SHAPE FUNCTIONS AND THEIR DERIVATIVES 

% AT THE GAUSS POINTS AND NODES. 

shape; 

o 

o 

% INPUT, PROCESS AND OUTPUT THE BOUNDARY CONDITIONS, 
bcs ; 

o, 

o 

% FORM THE COEFFICIENT MATRIX AND APPLY THE BOUNDARY CONDITIONS, 
frmtrx; 

o 

o 

% SOLVE THE LINEAR EQUATIONS. 
neqn=2*nnp; 

[uv, iflag] =elimin (a, neqn) ; 
if (iflag == 1) 

fprintf (fid6, ['\n', ... 

'MATRIX ILL-CONDITIONING DETECTED IN EQUATION SOLVER','\n']); 
error('MATRIX ILL-CONDITIONING DETECTED IN EQUATION SOLVER'); 
end; 

o, 

o 

% OUTPUT THE NODAL DISPLACEMENTS, ELEMENT STRESSES AND FORCES ON THE 
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% BOUNDARY SEGMENTS, 
output; 

o 

o 

end %program BEM2PQ 


function intitle(varargin) 

o 

o 

% SUBPROGRAM TO INPUT PROBLEM TITLE AND WHETHER PLANE 
% STRESS OR PLANE STRAIN. ALSO THE MATERIAL PROPERTIES. 

o, 

o 

global nu e fid5 fid6 

o, 

o 

% INPUT THE PROBLEM TITLE. 
title=fgetl (fid5) ; 

fprintf (fid 6, [ 'QUADRATIC BOUNDARY ELEMENT SOLUTION FOR', ... 

' TWO DIMENSIONAL ELASTIC PROBLEM\n\n%s\n'], title); 

o, 

o 

% INPUT THE PROBLEM CASE TYPE: 

% 'STRESS' DEFINES PLANE STRESS 

% 'STRAIN' DEFINES PLANE STRAIN 

% THE DEFAULT IS PLANE STRESS. 
pcase=fgetl (fid5) ; 

fprintf (fid 6 , [ ' \n ' , 'UNDER PLANE ' , ' %s ' , ' CONDITIONS ' , ' \n ' ] , pease) ; 

o_ 

o 

% INPUT AND OUTPUT YOUNG'S MODULUS AND POISSON'S RATIO. 
line=fgetl (fid5) ; 
vec=str2num(line); 
e=vec(l); nu=vec(2); 

fprintf (fid 6 , [ ' \n ' , ' YOUNGS MODULUS = ', , %12.4e', ... 

' POISSONS RATIO = ','%5.3f','\n'],e,nu); 

o 

o 

% MODIFY PROPERTIES IF CASE IS ONE OF PLANE STRAIN. 
strain='STRAIN'; 
if(strempi(pease,strain) == 1) 
e=e/(1.-nu A 2); 
nu=nu/(1.-nu); 
end; 

o 

o 
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return; 

end %function intitle 


function meshq(varargin) 

o, 

o 

% SUBPROGRAM TO READ IN AND GENERATE THE GEOMETRIC DATA FOR A MESH OF 
% QUADRATIC ELEMENTS. 

o 

o 

global node maxi ynode xnode nnp nel neend m3 ml angstore yeend xeend ... 
isegelem mfirst isegend ysend xsend mlast maxnel nelb nsegb . . . 
nsegtot nbound maxnb fid5 fid6 

o, 

o 

% INPUT THE NUMBER OF SEPARATE BOUNDARIES. 
line=fgetl (fid5) ; 
nbound=str2num(line) ; 

o, 

o 

% TEST THE NUMBER OF BOUNDARIES, 
if(nbound < 1 || nbound > maxnb) 

fprintf (fid6, ['in', ’ NBOUND = * , ’ %4i ’ , ... 

’ OUTSIDE PERMITTED RANGE 1 TO%4i\n'],nbound,maxnb); 
error(’nbound error’); 
end; 

o, 

o 

% FOR EACH BOUNDARY IN TURN INPUT THE NUMBER OF SEGMENTS. 

nel=0; 

ieend=0; 

nsegtot=0; 

mmaxb=0; 

for ibound=l:nbound; % each boundary in turn 
nelb(ibound)=0; 
mminb=mmaxb+l; 
line=fgetl (fid5) ; 
nsegb(ibound)=str2num(line); 
nsegtot=fix (nsegtot+nsegb (ibound) ) ; 

o, 

o 

% TEST THE NUMBER OF SEGMENTS. 

if(nsegtot < 1 || nsegtot > maxnel) 

fprintf (fid 6 , [ ’ \n ’ , ’ NSEGTOT =’,’%6i’, ... 

’ OUTSIDE PERMITTED RANGE 1 TO’,’%6i’,’\n '], nsegtot f maxnel); 
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error('nsegtot error'); 


end; 


o 

o 


% INPUT THE CARTESIAN GLOBAL COORDINATES OF THE END POINTS OF THE 
% SEGMENTS. TAKE THE END POINTS CONSECUTIVELY, KEEPING THE DOMAIN 
% TO THE LEFT OF THE DIRECTION OF NUMBERING. 

% Each line of input data is assumed to contain an arbitrary number 
% of pairs of x,y coordinates. 
isend=0; 

while isend < nsegb(ibound); 
line=fgetl (fid5) ; 
vec=str2num(line) ; 
veclength=length(vec); 
np=veclength/2 ; 
for ip=l:np; 
isend=isend+l; 
xsend(isend)=vec (2*ip-l); 
ysend (isend) =vec (2*ip) ; 
end; 
end; 
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o, 

o 

% DEFINE THE FIRST END POINT ON THE CURRENT BOUNDARY. 
ieend=ieend+l; 
xeend(ieend)=xsend(1); 
yeend(ieend)=ysend(1) ; 

o, 

o 

% FOR EACH OF THE SEGMENTS (BETWEEN ENDS 1 AND 2, 2 AND 3, ETC.) 

% INPUT THE RADIUS OF CURVATURE (+VE FOR CONVEX WITH CENTRE OF 
% CURVATURE INSIDE DOMAIN, -VE FOR CONCAVE), THE NUMBER OF 
% ELEMENTS IN THE SEGMENT, AND THE LENGTH RATIO BETWEEN SUCCESSIVE 
% ELEMENTS IN THE DIRECTION OF NUMBERING. 
isegmax=nsegtot; 

isegmin=isegmax-nsegb(ibound)+1; 

for iseg=isegmin:isegmax; % each segment in turn 
line=fgetl (fid5) ; 
vec=str2num(line) ; 

rseg=vec(l); nelseg=vec(2); ratseg=vec(3); 

o, 

o 

% FIND AND TEST THE NUMBER OF ELEMENTS SO FAR. 
nel=fix (nel+nelseg) ; 

nelb (ibound) =fix (nelb (ibound) +nelseg) ; 
if (nel < 1 | | nel > maxnel) 

fprintf (fid 6, ['\n', 'NEL =','%6i',' OUTSIDE PERMITTED' ... 

' RANGE 1 TO','%6i','\n'],nel,maxnel); 
error('nel error'); 
end; 

o 

o 

% FIRST AND LAST ELEMENTS ON CURRENT SEGMENT, 
mlast (iseg) =fix (nel) ; 
mfirst (iseg) =fix (nel-nelseg+1) ; 
mmaxb=mmaxb+nelseg; 

o, 

o 

% COORDINATES OF THE FIRST END POINT OF THE SEGMENT. 
isend=iseg-isegmin+l; 
xfirst=xsend (isend) ; 
yfirst=ysend (isend) ; 

o, 

o 

% COORDINATES OF THE LAST END POINT OF THE SEGMENT. 
isend=isend+l; 
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o, 

o 


o, 

o 


o 

o 


o 

o 


o 

o 


o, 

o 


if(iseg == isegmax) 
isend=l; 
end; 

xlast=xsend(isend); 
ylast=ysend(isend); 

GENERATE ELEMENT DATA FOR A STRAIGHT SEGMENT, 
if(rseg == 0.) 

DEFINE THE ELEMENT END POINT COORDINATES ON THE SEGMENT, 
for m=l:nelseg; % each element in turn 
ieend=ieend+l; 
isegend (ieend) =fix (iseg) ; 
if (ratseg == 1. ) 

xeend (ieend) =xfirst+ (xlast-xfirst) * (m) / (nelseg) ; 
yeend (ieend) =yfirst+ (ylast-yfirst) * (m) / (nelseg) ; 
end; 

if (ratseg ~= 1. ) 

xeend (ieend) =xfirst+ (xlast-xfirst) * (1. -ratseg A m) / . . . 
(1.-ratseg A nelseg); 

yeend (ieend) =yfirst+ (ylast-yfirst) * (1. -ratseg A m) / . . . 
(1.-ratseg A nelseg); 

end; 

end; % each element in turn 

DEFINE THE ELEMENT NODES AND COORDINATES. 
ieend=ieend-nelseg; 

for ielseg=l:nelseg; % each element in turn 
ieend=ieend+l; 
m=mfirst (iseg) +ielseg-l; 
il=2 *m-l; 
i2=il+l; 
i3=il+2; 

if(iseg == isegmax && ielseg == nelseg) 
i3=node(mminb,1); 
end; 

node (m, 1) =fix (il) ; 
node (m, 2) =fix (12) ; 
node (m, 3) =fix (13) ; 
isegelem (m) =fix (iseg) ; 


o, 

o 


o, 

o 


o 

o 


o 

o 


o 

o 


o, 

o 
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if(iseg == isegmin && ielseg == 1) 
xnode(il)=xeend(ieend-1); 
ynode(il)=yeend(ieend-1); 
end; 

xnode(13)=xeend(ieend) ; 

ynode(13)=yeend(ieend); 

xnode(i2)=0.5*(xnode(il)+xnode(i3)); 

ynode(i2)=0.5*(ynode(il)+ynode(i3) ) ; 

o 

o 

% STORE THE NUMBERS OF THE ADJACENT ELEMENTS, 
ml (m) =fix (m-1) ; 
m3 (m) =fix (m+1) ; 
end; % each element in turn 
end; 

o 

o 

% GENERATE ELEMENT DATA FOR A SEGMENT IN THE FORM OF A CIRCULAR ARC. 
if(rseg ~= 0.) 

o, 

o 

% LOCATE THE CENTRE OF THE ARC. 
xmid= (xfirst+xlast) /2 . ; 
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ymid= (yfirst+ylast) /2 . ; 

alseg=sqrt ( (xlast-xfirst) A 2+ (ylast-yfirst) A 2) ; 
alperp2=rseg A 2-(alseg/2.) A 2; 
if(abs (alperp2) < 1.0e-6*rseg A 2) 
alperp2=0.; 
end; 

if (alperp2 < -1.0e-6*rseg A 2) 

fprintf (fid6, [ ' \n ' , ' DATA ERROR FOR SEGMENT NUMBER6i 1 , ... 

' \n ' ,’NOT POSSIBLE TO CREATE A CIRCULAR ARC’,’\n'],iseg); 
error('circular arc error'); 
end; 

alperp=sqrt(alperp2); 
uvflx= (xlast-xfirst) /alseg; 
uvfly= (ylast-yfirst) /alseg; 
fact=l.; 
if(rseg < 0.) 

fact=-l.; 
end; 

xcent=xmid-alperp ,lr uvfly ,,r f act; 
ycent=ymid+alperp ,lr uvflx ,,r f act; 

o 

o 

% FIND THE ANGLE SUBTENDED THERE BY THE SEGMENT, 
if(alperp ~= 0.) 

angseg=2.*atan(alseg^O.5/alperp); 
end; 

if(alperp == 0.) 

angseg=pi; 
end; 

o 

o 

% DEFINE THE ELEMENT END POINT COORDINATES ON THE SEGMENT. 
angfir=atan2 (yfirst-ycent, xfirst-xcent) ; 
angstore(ieend)=0. ; 

for m=l:nelseg; % each element in turn 
ieend=ieend+l; 
isegend (ieend) =fix (iseg) ; 
if (ratseg == 1.) 

ang=angseg*(m)/(nelseg) ; 
end; 

if(ratseg ~= 1.) 

ang^angseg^(1.-ratseg A m)/(1.-ratseg A nelseg); 


150 


Download free eBooks at bookboon.com 


Boundary Element Methods for Engineers: 
Part 2: Plane Elastic Problems 


Appendix E: Matlab Version of Quadratic Boundary 
Element Program for Plane Elastic Problems 


end; 

if (rseg < 0.) 

ang=-ang; 

end; 

xeend (ieend) =xcent+abs (rseg) *cos (angfir+ang) ; 
yeend (ieend) =ycent+abs (rseg) *sin (angfir+ang) ; 
angstore(ieend)=ang; 
end; % each element in turn 

o 

o 

% DEFINE THE ELEMENT NODES AND COORDINATES. 
ieend=ieend-nelseg; 

for ielseg=l:nelseg; % each element in turn 
ieend=ieend+l; 
m=mfirst (iseg) +ielseg-l; 
il=2*m-l; 
i2=il+l; 
i3=il+2; 

if(iseg == isegmax && ielseg == nelseg) 
i3=node(mminb,1); 
end; 

node (m, 1) =fix (il) ; 
node (m, 2) =fix (i2) ; 
node (m, 3) =fix (i3) ; 
isegelem (m) =fix (iseg) ; 
if(iseg == isegmin && ielseg == 1) 
xnode(il)=xeend(ieend-1); 
ynode(il)=yeend(ieend-1); 
end; 

xnode(i3)=xeend(ieend) ; 
ynode(13)=yeend(ieend); 

ang=0.5^(angstore(ieend-1)+angstore(ieend)) ; 
xnode (12) =xcent+abs (rseg) *cos (angfir+ang) ; 
ynode (12) =ycent+abs (rseg) *sin (angfir+ang) ; 

o 

o 

% STORE THE NUMBERS OF THE ADJACENT ELEMENTS, 
ml (m) =fix (m-1) ; 
m3 (m) =fix (m+1) ; 
end; % each element in turn 
end; 

o 

o 
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end; % each segment in turn 

o 

o 

% ADJACENT ELEMENTS FOR END ELEMENTS OF THE BOUNDARY, 
ml (mminb) =fix (mmaxb) ; 
m3 (mmaxb) =fix (mminb) ; 
end; % each boundary in turn 
neend=fix (ieend) ; 
nnp=fix (nel*2) ; 

o 

o 

% DETERMINE THE MAXIMUM DIMENSION OF THE SOLUTION DOMAIN. 
maxl=0.; 

for i=l:nnp; % each node in turn 

for j=l:nnp; % each other node in turn 

dist=sqrt((xnode(i)-xnode(j)) A 2+(ynode(i)-ynode(j)) A 2) ; 
if(dist > maxi) 
maxl=dist; 
end; 

end; % each other node in turn 
end; % each node in turn 

o, 

o 
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return; 

end %function meshq 


function mshout(varargin) 

o, 

o 

% SUBPROGRAM TO WRITE OUT THE MESH DATA. 

o, 

o 

global node maxi ynode xnode nnp nel fid7 

o 

o 

% OUTPUT THE NUMBERS OF ELEMENTS AND NODES, ALSO THE NODAL 
% COORDINATES. 

fprintf (fid7, ['GEOMETRIC DATA FOR THE MESH ' , ' \n ' , ' \n ' , ... 

' NUMBER OF ELEMENTS =','%6i','\n','\n', ... 

' NUMBER OF NODAL POINTS =','%6i’,’\n','\n', .. 

'COORDINATES OF THE NODAL POINTS','\n','\n' , ... 

' I X Y ' , . . . 

' IX Y '],nel,nnp); 

for i=l:nnp; 

if(mod(i,2) == 1) ; 

fprintf (fid7, [ ' \n ' , ' %6i ' , ' %12.4e' , ' %12.4e' ] , ... 

i, xnode(i),ynode(i)); 

end; 

if(mod(i,2) == 0); 

fprintf (fid7, [ ' %6i ' , ' %12.4e ' , ' %12.4e ' ] , ... 

i, xnode(i),ynode(i)); 

end; 

end; 


% OUTPUT THE ELEMENT NODE NUMBERS. 

fprintf (fid7 , [ ' \n ' , ' \n ' , ' ELEMENT NODE NUMBERS ' , ' \n ' , ' \n ' , 


' M ND1 

ND2 

ND3', ... 

' M ND1 

ND2 

ND3']); 

for m=l:nel; 

if(mod(m, 2) == 1) 

fprintf (fid7, [ ' \n ' , ' 


' , ' %5i' ] , m) ; 

end; 

if(mod(m, 2) == 0) 

fprintf (fid7, [ ' 

! ! O, 

r ° 

5i'],m); 


end; 
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for in=l:3; 

fprintf (fid7, f %5i ’ , node (m, in) ) ; 
end; 
end; 

o 

o 

% SCALE THE NODAL POINT COORDINATES, 
for i=l:nnp; % each node in turn 
xnode(i)=xnode(i)/maxi; 
ynode(i)=ynode(i)/maxi; 
end; % each node in turn 

o 

o 

return; 

end %function mshout 


function bcs(varargin) 

o 

o 

% SUBPROGRAM TO INPUT, PROCESS AND OUTPUT THE BOUNDARY CONDITIONS. 

o, 

o 

global nodepc idirpc nnp nbcpc nbcu e estore node tmy ty tmx tx ... 
ibcn ibce maxi velem v uelem u m3 ml nel mlast mfirst sigsnseg . . . 
signnseg sigsn signn unmx unmy nsegtot isegbc nbcs vseg useg ... 
maxnpc maxnel nbct ungy ungx fid5 fid6 

o, 

o 

% FIRST FIND THE UNIT NORMAL COMPONENTS AT THE ELEMENT NODES, 
for m=l:nel; % each element in turn 

for in=l:3; % each element node in turn 
igauss=in+10; 
it=l; 
ic=l ; 

jacobi(m,igauss,it,ic); 
unmx(m,in)=ungx; 
unmy(m,in)=ungy; 

end; % each element node in turn 
end; % each element in turn 

o, 

o 

% INPUT THE NUMBERS OF SEGMENTS SUBJECT TO EACH TYPE OF BOUNDARY 
% CONDITION. ALSO THE NUMBER OF POINT CONSTRAINTS. 

% NBCU - PRESCRIBED DISPLACEMENTS. 

% NBCS - NON-ZERO PRESCRIBED STRESSES. 
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% ANY SEGMENT NOT INCLUDED IS ASSUMED TO BE SUBJECT TO ZERO 
% STRESSES. 

% NBCPC - NODAL POINT DISPLACEMENT CONSTRAINTS. 

o 

o 

line=fgetl (fid5) ; 
vec=str2num(line); 

nbcu=vec(1); nbcs=vec(2); nbcpc=vec(3); 

o, 

o 

% TEST THESE BOUNDARY CONDITION NUMBERS. 
nbct=fix (nbcu+nbcs) ; 

if(nbcu < 0 || nbcu > maxnel || nbcs < 0 || nbcs > maxnel ... 

|| nbct < 0 || nbct > maxnel) 

fprintf (fid6 , [ '\n' , 'NBCU =' , '%6i' , ' NBCS = ' , ' %6i ' , ' \n ' , ... 

' NBCT = ' ,'%6i',' OUTSIDE PERMITTED RANGE 0 TO', ... 

'%6i 'f '\n'],nbcu,nbcs , nbct,maxnel); 
error('numbers of boundary conditions error'); 
end; 

if(nbcpc < 0 || nbcpc > maxnpc) 

fprintf (fid6, ['\n', ' NBCPC = ' , ' %4i ' , ... 

' OUTSIDE PERMITTED RANGE 0 TO','%3i\n'],nbcpc,maxnpc); 
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error('number of point constraints error'); 
end; 

o 

o 

% INITIALISE THE BOUNDARY CONDITION STORAGE ARRAYS, 
for m=l:nel; % each element in turn 
ibce(m)=2; 

for in=l:3; % each element node in turn 
signn(m,in)=0.; 
sigsn(m,in)=0.; 
tmx(m,in)=0.; 
tmy(m, in)=0 . ; 

end; % each element node in turn 
end; % each element in turn 

o_ 

o 

% INPUT, STORE AND OUTPUT THE PRESCRIBED DISPLACEMENT CONDITIONS, 
if(nbcu > 0) 

for ibcu=l:nbcu; 
line=fgetl (fid5) ; 
vec=str2num(line) ; 

isegbc(ibcu)=vec(1); useg(ibcu)=vec(2); vseg(ibcu)=vec(3); 
end; 

fprintf (fid6, [ '\n' , 'PRESCRIBED DISPLACEMENT BOUNDARY CONDITIONS','\n']); 
for ibcu=l:nbcu; % each segment with prescribed displacements 
iseg=isegbc(ibcu); 
if(iseg < 1 || iseg > nsegtot) 

fprintf (fid6, [' \n' ,' ISEG = ','%6i',' OUTSIDE PERMITTED RANGE' ... 

' 1 TO','%6i','\n'],iseg, nsegtot) ; 
error('iseg error'); 
end; 

for m=mfirst (iseg) :mlast (iseg) ; % each element on current segment 
ibce(m)=1; 

uelem(m)=useg(ibcu); 
velem(m)=vseg (ibcu); 
end; % each element on current segment 

o 

o 

fprintf (fid6, ['\n', 'U =', '%12.4e', ' V =', '%12.4e', . . . 

' ON ELEMENTS','%4i', ' TO', '%4i', '\n'] , ... 

useg (ibcu) , vseg (ibcu) , mfirst (iseg) , mlast (iseg) ) ; 
end; % each segment with prescribed displacements; 
end; 

o 

o 
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% INPUT, STORE AND OUTPUT THE PRESCRIBED STRESS CONDITIONS, 
if(nbcs > 0) 

for ibcs=l:nbcs; 
line=fgetl (fid5) ; 
vec=str2num(line) ; 

isegbc(ibcs)=vec(1)/ signnseg(ibcs)=vec(2); sigsnseg(ibcs)=vec(3) ; 
end; 

fprintf (fid6, [ '\n' , 'PRESCRIBED STRESS BOUNDARY CONDITIONS','\n']); 
for ibcs=l:nbcs; % each segment with prescribed stresses 
iseg=isegbc(ibcs); 
if(iseg < 1 || iseg > nsegtot) 

fprintf (fid6, ['\n','ISEG = ','%6i',' OUTSIDE PERMITTED RANGE' ... 

' 1 TO%6i\n'],iseg,nsegtot); 
error('iseg error'); 
end; 

for m=mfirst (iseg) :mlast (iseg) ; % each element on current segment 
ibce(m)=2; 

o, 

o 

% FIND THE TRACTIONS AT THE ELEMENT NODES FROM THE PRESCRIBED STRESSES. 
% ALSO STORE THE STRESSES AT THE ELEMENT NODES, 
for in=l:3; % each element node in turn 

tmx(m,in)=signnseg(ibcs)*unmx(m,in)-sigsnseg(ibcs)*unmy(m,in); 
tmy(m,in)=signnseg(ibcs)*unmy(m,in)+sigsnseg(ibcs)*unmx(m,in); 
signn(m,in)=signnseg(ibcs); 
sigsn(m,in)=sigsnseg(ibcs); 
end; % each element node in turn 
end; % each element on current segment 

o 

o 

fprintf (fid6, ['\n', 'SIGNN = ' , ' %12.4e ' , ' SIGSN =','%12.4e', ... 

' ON ELEMENTS','%4i', ' TO', '% 4i', '\n'] , ... 

signnseg (ibcs) , sigsnseg (ibcs) ,mfirst (iseg) ,mlast (iseg) ) ; 
end; % each segment with prescribed stresses 
end; 

o 

o 

% ASSEMBLE THE VECTOR OF KNOWN VARIABLES, APPLYING SCALING. 

% FIRST INITIALISE THE NODAL DISPLACEMENTS AND TRACTIONS, 
for i=l:nnp; % each node in turn 
u (i) =0 .; 
v (i) =0 . ; 
tx(i)=0.; 
ty(i)=0.; 
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end; % each node in turn 

for m=l:nel; % each element in turn 

for in=l:3; % each element node in turn 
i=node(m,in); 
ibcn (i) =fix (ibce (m) ) ; 

o, 

o 

% CHECK THE CONDITION APPLIED TO THE ADJACENT ELEMENT FOR AN ELEMENT 
% end NODE. IF THE CONDITION IS PRESCRIBED DISPLACEMENTS BUT THE 
% EXISTING CONDITION AT THE NODE IS NOT, THEN IMPOSE PRESCRIBED 
% DISPLACEMENTS, 
madj =m; 
if(in == 1) 
madj=ml(m); 
end; 

if(in == 3) 
madj =m3(m); 
end; 

if (ibce (madj ) == 1 && ibcn(i) == 2) 
ibcn(i)=1; 
end; 
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o, 

o 

% STORE KNOWN VARIABLES, 
if (ibcn (i) == 1) 
if(ibce(m) == 1) 
if (u (i) == 0.) 

u(i)=uelem(m)/maxi; 
end; 

if(v(i) == 0.) 

v(i)=velem(m)/maxi; 
end; 

if(u(i) ~= 0.) 

u(i)=0.5*(u(i)+uelem(m)/maxi); 
end; 

if(v(i) — 0.) 

v(i)=0.5*(v(i)+velem(m)/maxi); 
end; 
end; 
end; 

if(ibce(m) == 2 && ibcn(i) == 2) 
if(tx(i) — 0.) 

tx(i)=tmx(m,in)/e; 
end; 

if(ty(i) == 0.) 

ty(i)=tmy(m,in)/e; 
end; 

if(tx(i) 0.) 

tx(i)=0.5*(tx(i)+tmx (m, in) /e ); 
end; 

if(ty(i) -= 0.) 

ty(i)=0.5*(ty(i)+tmy(m,in)/e); 
end; 
end; 

end; % each element node in turn 
end; % each element in turn 

o 

o 

% SCALE YOUNG'S MODULUS. 

estore=e; 

e=l.; 

o 

o 
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% INPUT, STORE AND OUTPUT NODAL POINT DISPLACEMENT CONSTRAINTS. 

% IDIRPC STORES DIRECTIONS OF THE CONSTRAINTS - 1 FOR X, 2 FOR Y. 

% SUCH CONSTRAINTS SHOULD NOT REQUIRE POINT FORCES TO MAINTAIN THEM, 
if(nbcu == 0 && nbcpc < 3) 

fprintf (fid6, [' \n' , ’ INSUFFICIENT DISPLACEMENT BOUNDARY CONDITIONS', ... 

' OR POINT CONSTRAINTS','\n']); 
error('displacement constraint error'); 
end; 

if(nbcpc > 0) 

for ibcpc=l:nbcpc; 
line=fgetl (fid5) ; 
vec=str2num(line) ; 

nodepc(ibcpc)=vec(1); idirpc(ibcpc)=vec(2); 
end; 

for ibcpc=l:nbcpc; % each point constraint in turn 
if(nodepc(ibcpc) < 0 || nodepc(ibcpc) > nnp) 

fprintf (fid6, ['\n', 'POINT CONSTRAINT NODE','%4i', ... 

' OUTSIDE RANGE 0 TO ','%4i','\n'],nodepc(ibcpc),nnp); 
error('point constraint node error'); 
end; 

if(idirpc(ibcpc) ~= 1 && idirpc(ibcpc) ~= 2) 

fprintf (fid6, ['\n', 'POINT CONSTRAINT DIRECTION','%3i ' , ... 

' FOR NODE','%4i',' OUTSIDE RANGE 1 TO 2','\n'], ... 

idirpc(ibcpc),nodepc(ibcpc)); 
error('point constraint direction error'); 
end; 

if(idirpc(ibcpc) == 1) 

fprintf (fid6, ['\n', 'NODE', '%4i', ' CONSTRAINED WITH U = O', ... 

'\n'],nodepc(ibcpc)); 

end; 

if(idirpc(ibcpc) == 2) 

fprintf (fid6, [ '\n' , 'NODE' , ' %4i ' , ' CONSTRAINED WITH V = O', ... 

'\n'],nodepc(ibcpc)); 

end; 

end; % each point constraint in turn 
end; 

o, 

o 

return; 

end %function bcs 
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function frmtrx(varargin) 

o, 

o 

% SUBPROGRAM TO FORM THE COEFFICIENT MATRIX AND RIGHT HAND SIDE VECTOR, 
% MODIFIED TO SUIT THE BOUNDARY CONDITIONS. 

o 

o 

global a idirpc nodepc nbcpc node v arowy u arowx ibcn estore tmy ... 
browy tmx browx ibce nel nnp 

o, 

o 

% DEFINE THE NUMBER OF COLUMNS IN THE EXTENDED COEFFICIENT MATRIX A. 
jmax=2 *nnp+l; 

o 

o 

% FORM THE COEFFICIENT MATRIX A, AND THE RIGHT HAND SIDE VECTOR B*T. 
for i=l:nnp; % take each node in turn as P 

o 

o 

% INITIALISE THE TWO ROWS OF MATRIX A CORRESPONDING TO THE NODE, 
for j=l:jmax; % each pair of coefficients in turn 
a(2*i-l,j)=0.; 
a(2*i,j)=0.; 

end; % each pair of coefficients in turn 

o, 

o 
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% INITIALISE THE ELEMENT CONTRIBUTIONS TO THE CURRENT ROWS OF THE 
% A AND B MATRICES. 

for m=l:nel; % each element in turn 

for in=l:3; % each element node in turn 
arowx(m,2*in-l)=0.; 
arowx(m,2*in)=0.; 
arowy(m,2*in-l)=0.; 
arowy(m,2*in)=0.; 
browx(m,2*in-l)=0.; 
browx(m,2*in)=0.; 
browy(m,2*in-l)=0.; 
browy(m,2*in)=0.; 
end; % each element node in turn 
end; % each element in turn 


o 

o 


% SET UP THE LOOP TO INTEGRATE OVER EACH ELEMENT IN TURN, 
for m=l:nel; % each element in turn 


o 

o 


% INTEGRATE THE KERNEL PRODUCTS OVER THE CURRENT ELEMENT, 
intker(i,m); 

end; % each element in turn 


o 

o 


% EVALUATE THE COEFFICIENTS AT THE DIAGONAL OF MATRIX A. 
aiixx=0.; 
aiixy=0.; 
aiiyx=0.; 
aiiyy=0.; 

for m=l:nel; % each element in turn 

for in=l:3; % each element node in turn 
aiixx=aiixx-arowx(m f 2*in-l); 
aiixy=aiixy-arowx(m,P^in); 
aiiyx=aiiyx-arowy(m f 2*in-l); 
aiiyy=aiiyy-arowy(m,2*in); 
end; % each element node in turn 
end; % each element in turn 
if(ibcn (i) == 2) 

a(2*i-l,2*i-l)=aiixx; 
a(2*i-l,2*i)=aiixy; 
a(2^i f 2^i-l)=aiiyx; 
a(2 *i , 2 *i)=aiiyy; 
end; 
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o, 

o 

% INITIALISE THE B*T VECTOR COEFFICIENTS. 
btx=0.; 
bty=0.; 

if (ibcn (i) == 1) 

btx=-aiixx*u(i)-aiixy*v(i); 
bty=-aiiyx*u (i) -aiiyy*v (i) ; 
end; 

o 

o 

% APPLY THE BOUNDARY CONDITIONS TO THE CURRENT ROWS OF A AND B, BY 
% CONSIDERING EACH ELEMENT IN TURN, 
for m=l:nel; % each element in turn 

for in=l:3; % each element node in turn 
j =node(m,in) ; 

o 

o 

% IF DISPLACEMENTS ARE PRESCRIBED OVER THE ELEMENT, 

% INTERCHANGE THE A AND B COEFFICIENTS, 
if(ibce(m) == 1) 

a (2*i-l,2*j-1)=a(2*i-l,2*j-1)-browx(m,2*in-l); 
a(2*i-l,2*j)=a (2*i-l,2*j)-browx(m,2*in); 
a(2 *i,2*j-1)=a(2 *i,2 *j-1) -browy (m, 2 *in-l) ; 
a(2*i,2*j)=a(2*i,2*j) -browy (m,2*in) ; 
btx=btx-arowx(m,2 *in-l)*u(j)-arowx(m,2 *in)*v(j); 
bty=bty-arowy(m,2^in-1)^u(j)-arowy(m,2^in)*v(j); 
end; 

o, 

o 

% IF STRESSES ARE PRESCRIBED OVER THE ELEMENT, STORE THE A MATRIX 
% COEFFICIENTS, EXCEPT AT A CORNER NODE WHERE PRESCRIBED DISPLACEMENTS 
% HAVE BEEN IMPOSED. 

if(ibce(m) == 2) 

btx=btx+(browx(m,2*in-l)*tmx(m,in)+browx(m,2*in)*tmy(m,in))/ 

estore; 

bty=bty+(browy(m,2*in-l)*tmx(m,in)+browy(m,2*in)^tmy(m,in))/ 

estore; 

if (ibcn (j) == 2) 

a (2*i-l,2*j-1)=a(2*i-l,2^j-1)+arowx(m,2*in-l); 
a(2^i-l,2 ,lr j)=a(2 ,,r i-l,2 ,,r j) +arowx (m, 2 *in) ; 
a(2*i,2*j-l)=a(2*i,2*j-l)+arowy(m,2 *in-l); 
a(2^i,2'* r j)=a(2 ,,r i,2 ,,r j) +arowy (m, 2 *in) ; 
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end; 

if (ibcn(j) == 1) 

btx=btx-arowx(m,2*in-l)*u(j)-arowx(m,2 *in)*v(j); 
bty=bty-arowy(m,2*in-l)*u(j)-arowy(m,2 *in)*v (j ); 
end; 
end; 

o 

o 

end; % each element node in turn 
end; % each element in turn 

o 

o 

% STORE THE B*T COEFFICIENTS AS EXTENSIONS OF MATRIX A. 
a(2*i-l,jmax)=btx; 
a (2*i,jmax)=bty; 

end; % take each node in turn as P 

o 

o 

% APPLY NODAL POINT DISPLACEMENT CONSTRAINTS, 
if (nbcpc > 0) 

for ibcpc=l:nbcpc; % each point constraint in turn 
irow=2*(nodepc(ibcpc)-1)+idirpc(ibcpc); 
for j=l:jmax; % each column in turn 
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a(irow,j)=0.; 

end; % each column in turn 
a(irow,irow) =1.; 

end; % each point constraint in turn 
end; 

o, 

o 

return; 

end %function frmtrx 


function shape(varargin) 

o, 

o 

% SUBPROGRAM TO EVALUATE AND STORE VALUES OF THE SHAPE FUNCTIONS AND 
% THEIR DERIVATIVES, AT THE GAUSS POINTS AND NODES. 

o 

o 

global sd sdl sfl egl ngauss sf zg wgl wg 

o, 

o 

% STORE APPROPRIATE COORDINATES AND WEIGHT FACTORS FOR NORMAL GAUSSIAN 
% QUADRATURE IN ARRAYS XG AND CG, ALSO THOSE FOR LOGARITHMIC FUNCTION 
% INTEGRATION IN XGL AND CGL. 
ngauss=8; 

zg(l)=-0.9602898564; 
zg(2)=-0.7966664774; 
zg(3)=-0.5255324099; 
zg(4)=-0.1834346424; 
zg(5)=-zg(4); 
zg(6)=-zg (3); 
zg(7)=-zg (2) ; 
zg(8)=-zg(1); 
wg(l)=0.1012285362; 
wg(2)=0.2223810344; 
wg (3)=0.3137066458; 
wg(4)=0.3626837833; 
wg(5)=wg(4) ; 
wg(6)=wg(3) ; 
wg(7)=wg (2) ; 
wg(8)=wg(1) ; 
egl(1)=0.013320244; 
egl(2)=0.079750429; 
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egl(3)=0.197871029; 
egl (4)=0.354153994; 
egl(5)=0.529458575; 
egl(6)=0.701814530; 
egl(7)=0.849379320; 
egl(8)=0.953326450; 
wgl(1)=0.164416605; 
wgl(2)=0.237525610; 
wgl(3)=0.226841984; 
wgl(4)=0.175754079; 
wgl(5)=0.112924030; 
wgl(6)=0.057872211; 
wgl(7)=0.020979074; 
wgl(8)=0.003686407; 

o 

o 

% NORMAL GAUSSIAN QUADRATURE. 

for igauss=l:ngauss; % each gauss point in turn 
zeta=zg(igauss) ; 

sf(1 , igauss)=0.5*zeta*(zeta-1. ) ; 
sf(2,igauss)= 1.-zeta A 2 ; 
sf (3, igauss)=0.5*zeta* (zeta+1.); 
sd(1 , igauss)=zeta-0.5; 
sd(2 , igauss)=-2.*zeta; 
sd(3,igauss)=zeta+0.5; 
end; % each gauss point in turn 

o, 

o 

% FOUR CASES OF LOGARITHMIC GAUSSIAN QUADRATURE TO CONSIDER. 


\—1 

II 

u 

H 
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TO 
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TO 
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TO 
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IC=4 
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OVER 

WHOLE 

ELEMENT 

FROM 

THIRD 

TO 

FIRST 

NODE 


o 

o 

for ic=l:4; % each case in turn 

for igauss=l:ngauss; % each gauss point in turn 
eta=egl(igauss) ; 
if(ic — 1) 

zeta=2.*eta-l.; 
end; 

if(ic == 2) 
zeta=eta; 
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end; 

if(ic == 3) 
zeta=-eta; 
end; 

if(ic == 4) 

zeta=l.-2.*eta; 
end; 

sfl (ic, 1, igauss) =0.5*zeta* (zeta-1. ) ; 
sfl (ic, 2, igauss) =l.-zeta A 2; 
sfl (ic, 3, igauss) =0.5*zeta* (zeta+1.) ; 
sdl(ic,1,igauss)=zeta-0.5; 
sdl(ic,2,igauss)=-2.*zeta; 
sdl(ic,3,igauss)=zeta+0.5; 
end; % each gauss point in turn 
end; % each case in turn 

o 

o 

% SHAPE FUNCTION DERIVATIVES AT THE NODES, STORED AS THOUGH THEY 
% ARE FOR GAUSS POINTS NUMBERED 11, 12 AND 13. 
for igauss=ll:13; % each element node in turn 
if(igauss == 11) 
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zeta=-l.; 
end; 

if(igauss == 12) 
zeta=0.; 
end; 

if(igauss == 13) 
zeta=l.; 
end; 

sd(1 , igauss)=zeta-0.5; 
sd (2 , igauss) =-2.*zeta; 
sd(3,igauss)=zeta+0.5; 
end; % each element node in turn 

o, 

o 

return; 

end %function shape 


function intker(i,m) 

o, 

o 

% SUBPROGRAM TO INTEGRATE THE KERNEL PRODUCTS FOR A PARTICULAR FORCE 
% POINT P (INDICATED BY NODE NUMBER I) OVER A PARTICULAR ELEMENT 
% (INDICATED BY ARGUMENT M) BY GAUSSIAN QUADRATURE. 

o, 

o 

global node jacob wgl browy browx sfl ngauss bk2yy wg bk2yx bk2xy ... 
bk2xx sf bkyy bkyx bkxy bkxx akyy arowy akyx akxy arowx akxx ... 
ungy ungx ynode xnode e nu 

o 

o 

% CONSTANT IN DISPLACEMENT KERNELS MULTIPLYING THE LOGARITHMIC TERM. 
cl=(1.+nu)*(3.-nu)/(4.*pi*e); 

o, 

o 

% COORDINATES OF POINT P. 
xp=xnode(i); 
yp=ynode(i); 

o 

o 

% SET UP THE ELEMENT NODE LOOP, 
for in=l:3; % each element node in turn 
j =node(m,in); 

o 

o 
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% IF P IS NOT THE CURRENT ELEMENT NODE, USE NORMAL GAUSSIAN QUADRATURE. 

if(i — j) 

for igauss=l:ngauss; % each gauss point in turn 

o 

o 

% EVALUATE JACOBIAN AND UNIT NORMAL VECTOR COMPONENTS, ALSO THE KERNELS 
% AT THE PARTICULAR GAUSS POINT FOR NORMAL QUADRATURE OVER THE WHOLE 
% ELEMENT. 

it=l ; 
ic=l; 

jacobi(m,igauss,it,ic); 

kernel(xp,yp,m,igauss,ungx,ungy); 

o, 

o 

% ACCUMULATE THE INTEGRALS. 
sfn=sf(in,igauss); 

arowx(m,2*in-l)=arowx(m,2*in-l)+wg(igauss)*akxx*sfn*jacob; 
arowx(m,2*in)=arowx(m,2* in)+wg(igauss)*akxy*sfn*j acob; 
arowy(m,2*in-l)=arowy(m,2*in-l)+wg(igauss)*akyx*sfn*jacob; 
arowy(m,2*in)=arowy(m,2*in)+wg(igauss)^akyy^sfn^j acob; 
browx(m,2*in-l)=browx(m,2*in-l)+wg(igauss)^bkxx^sfn^jacob; 
browx(m,2*in)=browx(m,2*in)+wg(igauss)*bkxy*sfn*j acob; 
browy(m,2*in-l)=browy(m,2*in-l)+wg(igauss)^bkyx^sfn^jacob; 
browy(m,2*in)=browy(m,2*in)+wg(igauss)^bkyy^sfn^j acob; 
end; % each gauss point in turn 
end; 

o, 

o 

% IF P IS THE CURRENT ELEMENT NODE, SOME LOGARITHMIC QUADRATURE 
% IS REQUIRED. 

if(i == j) 

o 

o 

% P AT THE FIRST NODE OF THE ELEMENT, 
if(in == 1) 

o, 

o 

% TERMS INVOLVING NORMAL QUADRATURE. 

for igauss=l:ngauss; % each gauss point in turn 
it=l; 
ic=l ; 

jacobi(m,igauss,it,ic); 
kern2(m,in,igauss); 
sfn=sf(in,igauss); 

browx(m,2*in-l)=browx(m,2*in-l)+wg(igauss)*bk2xx*sfn*jacob; 
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browx(m,2*in)=browx(m,2*in)+wg(igauss)*bk2xy*sfn*j acob; 
browy (m,2*in-l) =browy (m, 2*in-l) +wg (igauss) *bk2yx*sfn* j acob; 
browy(m,2*in)=browy(m,2*in)+wg(igauss)*bk2yy*sfn*j acob; 
end; % each gauss point in turn 

o 

o 

% TERMS INVOLVING LOGARITHMIC QUADRATURE. 

for igauss=l:ngauss; % each gauss point in turn 
it=2 ; 
ic=l; 

jacobi(m,igauss,it, ic); 
sfn=sfl (ic, in, igauss) ; 
dzde=2.; 

browx (m, 2*in-l) =browx (m, 2 ,lr in-l) +wgl (igauss) *cl*sfn^jacob^dzde; 
browy(m,2 *in)=browy(m,2 *in)+wgl(igauss)*cl*sf n *j acob*dzde; 
end; % each gauss point in turn 
end; 

o 

o 

% P AT THE SECOND NODE OF THE ELEMENT, 
if(in == 2) 

o, 

o 
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% TERMS INVOLVING NORMAL QUADRATURE. 

for igauss=l:ngauss; % each gauss point in turn 
it=l; 
ic=l; 

jacobi(m,igauss,it,ic); 
kern2(m,in,igauss); 
sfn=sf(in,igauss); 

browx(m,2*in-l)=browx(m,2*in-l)+wg(igauss)*bk2xx*sfn*jacob; 
browx(m,2*in)=browx(m,2* in)+wg(igauss)*bk2xy*sfn*j acob; 
browy(m,2*in-l)=browy(m,2*in-l)+wg(igauss)*bk2yx*sfn*jacob; 
browy (m, 2*in) =browy (m, 2* in) +wg (igauss) *bk2yy*sfn* j acob; 
end; % each gauss point in turn 

o, 

o 

% TERMS INVOLVING LOGARITHMIC QUADRATURE. 

for igauss=l:ngauss; % each gauss point in turn 
it=2 ; 
ic=2 ; 

jacobi(m,igauss,it,ic); 
sfn=sfl (ic, in, igauss) ; 
dzde=l.; 

browx(m,2*in-l)=browx(m,2*in-l)+wgl(igauss)*cl*sfn^jacob^dzde; 
browy(m,2 *in)=browy(m,2 *in)+wgl(igauss)*cl*sfn^j acob*dzde; 
ic=3 ; 

jacobi(m,igauss,it,ic); 
sfn=sfl (ic, in, igauss) ; 
dzde=l.; 

browx (m, 2 ,lr in-l) =browx (m, 2 ,lr in-l) +wgl (igauss) *cl*sfn^jacob^dzde; 
browy(m,2 *in)=browy(m,2 *in)+wgl(igauss)*cl*sfn^j acob*dzde; 
end; % each gauss point in turn 
end; 

o, 

o 

% P AT THE THIRD NODE OF THE ELEMENT, 
if (in == 3) 

o 

o 

% TERMS INVOLVING NORMAL QUADRATURE. 

for igauss=l:ngauss; % each gauss point in turn 
it=l ; 
ic=l ; 

jacobi(m,igauss,it,ic); 
kern2(m,in,igauss); 
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sfn=sf(in,igauss); 

browx (m, 2*in-1) =browx (m, 2*in-l) +wg (igauss) *bk2xx*sfn*jacob; 
browx(m,2* in)=browx(m,2* in)+wg(igauss)*bk2xy*sfn*j acob; 
browy(m,2*in-l)=browy(m,2*in-l)+wg(igauss)*bk2yx*sfn*jacob; 
browy(m,2*in)=browy(m,2* in)+wg(igauss)*bk2yy*sfn*j acob; 
end; % each gauss point in turn 

o 

o 

% TERMS INVOLVING LOGARITHMIC QUADRATURE. 

for igauss=l:ngauss; % each gauss point in turn 
it=2 ; 
ic=4 ; 

jacobi(m,igauss,it,ic); 
sfn=sfl (ic, in, igauss) ; 
dzde=2.; 

browx (m, 2*in-l) =browx (m, 2 ,lr in-l) +wgl (igauss) *cl*sfn^jacob^dzde; 
browy(m,2 *in)=browy(m,2 *in)+wgl(igauss)^cl^sfn^j acob*dzde; 
end; % each gauss point in turn; 
end; 
end; 

end; % each element node in turn 

o 

o 

return; 

end %function intker 


function jacobi(m,igauss,it,ic) 

o 

o 

% SUBPROGRAM TO EVALUATE THE JACOBIAN AND THE COMPONENTS OF THE UNIT 
% NORMAL VECTOR AT A PARTICULAR GAUSS POINT. 

% M INDICATES THE ELEMENT NUMBER. 

% IGAUSS INDICATES THE GAUSS POINT NUMBER. 

% IT INDICATES THE TYPE OF QUADRATURE. 

% IT=1 - NORMAL GAUSSIAN QUADRATURE. 

% IT=2 - LOGARITHMIC GAUSSIAN QUADRATURE. 

% IC INDICATES THE CASE NUMBER FOR LOGARITHMIC GAUSSIAN QUADRATURE, 

% AS DEFINED IN SUBROUTINE SHAPE. 

o 

o 

global jacob ungy ungx node ynode xnode sdl sd 
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o, 

o 

% CALCULATE THE DERIVATIVES OF THE GLOBAL COORDINATES WITH RESPECT TO 
% THE LOCAL INTRINSIC COORDINATE. 
dx=0.; 
dy=0.; 

for in=l:3; % each element node in turn 
if(it — 1) 

sfnd=sd(in,igauss); 
end; 

if(it == 2) 

sfnd=sdl(ic,in,igauss); 
end; 

j =node(m,in); 
dx=dx+sfnd*xnode(j ) ; 
dy=dy+sfnd*ynode(j); 
end; % each element node in turn 

o 

o 

% COMPONENTS OF THE LOCAL OUTWARD NORMAL VECTOR AT THE GAUSS POINT. 

ungx=dy; 

ungy=-dx; 
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o, 

o 

% JACOBIAN OF THE COORDINATE TRANSFORMATION. 
jacob=sqrt(ungx A 2+ungy A 2); 

o 

o 

% SCALE THE VECTOR COMPONENTS TO GIVE THE UNIT NORMAL VECTOR. 
ungx=ungx/j acob; 
ungy=ungy/j acob; 

o, 

o 

return; 

end %function jacobi 


function kernel(xp,yp, m,igauss,unx,uny) 

o 

o 

% SUBPROGRAM TO EVALUATE THE KERNELS WHEN P IS NOT THE CURRENT 
% ELEMENT NODE. 

% XP, YP INDICATE THE GLOBAL COORDINATES OF POINT P. 

% M INDICATES THE ELEMENT NUMBER. 

% IGAUSS INDICATES THE NUMBER OF THE GAUSS POINT, Q. 

o 

o 

global bkyy bkxy bkyx bkxx akyy akyx akxy akxx nu e node ynode xnode sf 

o 

o 

% COORDINATES OF GAUSS POINT Q. 
xq=0.; 
yq=0.; 

for in=l:3; % each element node in turn 
sfn=sf(in,igauss); 
j =node(m,in); 
xq=xq+sfn*xnode(j); 
yq = yq+sf n *ynode(j); 
end; % each element node in turn 

o 

o 

% COMPONENTS AND MAGNITUDE OF THE RADIUS VECTOR FROM P TO Q. 

rx=xq-xp; 

ry=yq-yp; 

rsq=rx A 2+ry A 2; 

r=sqrt(rsq); 

rx=rx/r; 

ry=ry/r; 
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o, 

o 

% RATE OF CHANGE OF R WITH THE NORMAL TO THE BOUNDARY AT Q. 
drdn=rx*unx+ry*uny; 

o 

o 

% PARAMETERS IN THE KERNEL FUNCTIONS. 

cl=-l./(4.*pi*r); 

c2=l.-nu; 

c3=2.*(1.+nu); 

c4=(1.+nu) A 2/(4.*pi*e); 

c5=(3.-nu)*log(1./r)/ (1.+nu); 

o 

o 

% EVALUATE THE KERNELS. 
akxx=cl*(c2+c3*rx*rx)*drdn; 
terml=c3*rx*ry*drdn; 
term2=c2* (rx*uny-ry*unx) ; 
akxy=cl*(terml-term2); 
akyx^cl^(terml+term2); 
akyy=cl*(02+03^ry^ry)^drdn; 
bkxx=c4 ,lr (c5 + rx*rx) ; 
bkxy=c4 *rx*ry; 
bkyx=bkxy; 
bkyy=c4*(c5+ry^ry); 

o 

o 

return; 

end %function kernel 


function kern2(m,in f igauss) 

o, 

o 

% SUBPROGRAM TO EVALUATE THE NON-SINGULAR LOGARITHMIC TERM IN THE 
% SECOND KERNEL WHEN P IS THE CURRENT ELEMENT NODE. 

% M INDICATES THE ELEMENT NUMBER. 

% IN INDICATES THE NUMBER OF THE ELEMENT NODE FORMING P. 

% IGAUSS INDICATES THE GAUSS POINT NUMBER. 

o, 

o 

global bk2yy bk2xy bk2yx bk2xx nu e ynode xnode zg node sf 

o, 

o 


% COORDINATES OF GAUSS POINT Q. 
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xq=0.; 
yq=0.; 

for ic=l:3; % each element node in turn 
sfn=sf(ic,igauss); 
j =node(m,ic); 
xq=xq+sfn*xnode(j); 
yq^yq+sfn^ynode(j); 
end; % each element node in turn 

o 

o 

% ELEMENT NODE NUMBERS. 
il=node(m,1); 
i2=node(m,2); 
i3=node(m,3); 

o 

o 

% EVALUATE THE INTRINSIC COORDINATE. 
zeta=zg(igauss); 

o 

o 

% P AT THE FIRST NODE OF THE ELEMENT, 
if(in == 1) 

xcomp=(zeta-2.)*xnode(il)+2.*(l.-zeta)*xnode(i2)+zeta*xnode(13); 
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ycomp=(zeta-2.)*ynode(il)+2.*(1.-zeta)*ynode(i2)+zeta*ynode (i3); 
rfn=sqrt(xcomp A 2+ycomp A 2); 

i=il; 

end; 

o 

o 

% P AT THE SECOND NODE OF THE ELEMENT, 
if(in == 2) 

xcomp=-0.5*(zeta-1.)*xnode(il)+zeta*xnode(i2)-0.5*(zeta+1.)*xnode(13) ; 
ycomp=-0.5*(zeta-1.)*ynode(il)+zeta*ynode(i2)-0.5*(zeta+1.)*ynode(i3); 
rfn=sqrt(xcomp A 2+ycomp A 2); 
i=i2; 
end; 

o, 

o 

% P AT THE THIRD NODE OF THE ELEMENT, 
if(in == 3) 

xcomp=-zeta*xnode(il)+2.*(l.+zeta)*xnode(i2)-(zeta+2.)*xnode(i3); 
ycomp=-zeta*ynode(il)+2.*(l.+zeta)*ynode(i2)-(zeta+2.)*ynode(i3) ; 
rfn=sqrt(xcomp A 2+ycomp A 2); 
i=i3; 
end; 

o 

o 

% COMPONENTS AND MAGNITUDE OF THE RADIUS VECTOR FROM P TO Q. 

xp=xnode(i); 

yp=ynode(i); 

rx=xq-xp; 

ry=yq-yp; 

rsq=rx A 2+ry A 2; 

r=sqrt(rsq); 

rx=rx/r; 

ry=ry/r; 

o, 

o 

% PARAMETERS IN THE KERNEL FUNCTIONS. 

c4=(1.+nu) A 2/(4.*pi*e); 

c5=(3.-nu)*log(1./rfn)/(1.+nu); 

o 

o 

% EVALUATE THE KERNELS. 
bk2xx=c4*(c5+rx*rx); 
bk2xy=c4 *rx*ry; 
bk2yx=bk2xy; 
bk2yy=c4*(c5 + ry*ry) ; 

o 

o 
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return; 

end %function kern2 


function output(varargin) 

o, 

o 

% SUBPROGRAM TO WRITE OUT THE NODAL POINT VALUES OF DISPLACEMENTS 
% AND ELEMENT STRESSES AND COMPUTE THE FORCES ON THE BOUNDARY SEGMENTS. 

o 

o 

global nsegtot fyseg fxseg node maxi tmy jacob wg tmx sf ngauss ... 

isegelem nel sige sigsn sigss signn nu e unmx unmy v sd u ty tx ... 
ibce estore nnp uv ibcn fid6 

o 

o 

% ARRANGE FOR U, V AND TX, TY TO CONTAIN THE NODAL DISPLACEMENTS 
% AND TRACTIONS, RESPECTIVELY, 
for i=l:nnp; % each node in turn 
if(ibcn (i) == 1) 
tx(i)=uv(2*i-l); 
ty(i)=uv(2*i); 
end; 

if (ibcn (i) == 2 ) 
u(i)=uv(2*i-l); 
v(i)=uv(2 *i); 
end; 

end; % each node in turn 

o 

o 


% HEADING FOR OUTPUT OF 

fprintf (fid6, [ 1 \n 1 , ’ NODAL 

NODAL DISPLACEMENTS 

POINT DISPLACEMENTS 

AND 

AND 

TRACTIONS. 

TRACTIONS' 

, '\n' 

,'\n'. 

' I U 

V 

TX 


TY' , 

'\n']); 


e=estore; 

for i=l:nnp; % each node in turn 

o 

o 

% REMOVE THE SCALING, 
u(i)=u(i)*maxl; 
v(i)=v(i)*maxl; 
tx(i)=tx(i)^estore; 
ty(i)=ty(i)Restore; 

o 

o 
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% OUTPUT THE NODAL VALUES OF DISPLACEMENTS AND TRACTIONS. 


fprintf (fid6, [ ' %4i 1 , 1 %15.6e f , ' %15.6e 1 , f %15.6e ’ , f %15.6e 1 , 
' \n ' ] , i, u (i) ,v(i) , tx (i) ,ty(i) ) ; 
end; % each node in turn 

o 

o 

% HEADING FOR OUTPUT OF ELEMENT STRESSES. 

fprintf (fid6, [ ' \n' , 'STRESSES AT THE NODES OF THE ELEMENTS' 
'\n','\n',' M IN SIGNN SIGSS 

' SIGE','\n 1 ]); 


SIGSN', ... 


o 

o 

% NORMAL AND SHEAR STRESSES AT THE NODES OF THE ELEMENTS, 
for m=l:nel; % each element in turn 

for in—1:3; % each element node in turn 
if(ibce(m) ~= 2) 

j =node(m,in); 
tmx(m,in)=tx (j) ; 
tmy(m,in)=ty(j); 

signn(m,in)=tmx(m,in)*unmx(m,in)+tmy(m,in)*unmy(m,in); 
sigsn(m,in)=-tmx (m, in)*unmy(m,in)+tmy(m,in)*unmx {m, in); 
end; 
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o, 

o 

% DIRECT STRESS ALONG THE BOUNDARY SURFACE. 
icl=node(m,1); 
ic2=node(m,2); 
ic3=node(m,3) ; 
igauss=in+10; 

dudz=sd(1,igauss)*u(icl)+sd(2,igauss)*u(ic2)+sd(3,igauss)*u (ic3); 
dvdz=sd (1, igauss) *v(icl)+sd(2, igauss) *v (ic2) +sd (3, igauss) *v (ic3) ; 
it=l; 
ic=l; 

jacobi(m,10+in,it,ic)/ 

ess=(-dudz*unmy(m,in)+dvdz*unmx(m,in))/(j acob*maxl); 
sigss(m,in)=e*ess+nu*signn(m,in); 

o, 

o 

% VON MISES EQUIVALENT STRESS. 


sige(m,in)=sqrt(signn(m,in) A 2+sigss(m,in) A 2 ... 
-signn(m,in)^sigss(m,in)+3.^sigsn(m,in) A 2); 


o 

o 


% OUTPUT THE STRESSES AT THE NODES OF THE ELEMENTS. 


fprintf (fid6, [ * %4i 1 , ' %4i 1 , ' %15.6e 1 , f %15.6e 1 , f %15.6e 1 , ’ %15.6e 1 , . . . 

’\n']in,signn(m,in),sigss(m,in),sigsn(m,in),sige(m,in)); 
end; % each element node in turn 
end; % each element in turn 

o, 

o 

% COMPUTE THE FORCES ON THE BOUNDARY SEGMENTS, 
for iseg=l:nsegtot; % each segment in turn 
fxseg(iseg)=0.; 
fyseg(iseg)=0.; 
end; % each segment in turn 
for m=l:nel; % each element in turn 
iseg=isegelem(m); 


o 

o 

% APPLY GAUSSIAN QUADRATURE. 
fxelem=0.; 
fyelem=0.; 

for in—1:3; % each element node in turn 

for igauss=l:ngauss; % each gauss point in turn 
sfn=sf(in,igauss); 
it=l; 
ic=l; 
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jacobi(m,igauss,it,ic); 

fxelem=fxelem+wg(igauss)*sfn*jacob*tmx(m,in)*maxl; 
fyelem=fyelem+wg(igauss)*sfn*j acob*tmy(m,in)*maxl; 
end; % each gauss point in turn 
end; % each element node in turn 

o, 

o 

% ACCUMULATE THE FORCES ON THE SEGMENT, 
fxseg(iseg)=fxseg(iseg)+fxelem; 
fyseg(iseg)=fyseg(iseg)+fyelem; 
end; % each element in turn 

o 

o 

% OUTPUT THE SEGMENT FORCE RESULTS. 

fprintf (fid6, [’\n', 'FORCES ACTING ON THE BOUNDARY SEGMENTS', ... 

'\n',’\n','SEGMENT FX FY','\n']); 

for iseg=l:nsegtot; 

fprintf (fid6, [ '%5i' , ’ %14.5e ' , ' %14.5e ' , ’ \n ' ] , ... 

iseg,fxseg(iseg),fyseg(iseg)); 

end; 

o. 

o 

return; 

end %function output 


function [x, iflag] =elimin (a, neqn) 

o_ 

o 

% SUBROUTINE FOR SOLVING SIMULTANEOUS LINEAR EQUATIONS BY GAUSSIAN 
% ELIMINATION WITH PARTIAL PIVOTING. 

o 

o 

% INITIALIZE ILL-CONDITIONING FLAG. 
iflag=0 ; 

for i=l:neqn; x(i)=0.; end; 

o 

o 

% SCALE EACH EQUATION TO HAVE A MAXIMUM COEFFICIENT MAGNITUDE OF UNITY. 
jmax=neqn+l; 

for i=l:neqn; % each equation in turn 
amax=0.; 

for j=l:neqn; % search for maximum 
absa=abs(a (i,j)); 


181 


Download free eBooks at bookboon.com 


Boundary Element Methods for Engineers: 
Part 2: Plane Elastic Problems 


Appendix E: Matlab Version of Quadratic Boundary 
Element Program for Plane Elastic Problems 


if(absa > amax) 
amax=absa; 
end; 

end; % search for maximum 

o 

o 

for j=l: jmax; % scale coefficients 
a(i,j)=a(i,j)/amax; 
end; % scale coefficients 

o 

o 

end; % each equation in turn 

o 

o 

% COMMENCE ELIMINATION PROCESS. 

for k=l:neqn-l; % eliminate each variable in turn 

o 

o 

% SEARCH LEADING COLUMN OF THE COEFFICIENT MATRIX FROM THE DIAGONAL 
% DOWNWARDS FOR THE LARGEST VALUE. 
imax=k; 

for i=k+l:neqn; % search for largest value 
if(abs (a (i, k) ) > abs(a(imax,k))) 
imax=i; 
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end; 

end; % search for largest value 

o 

o 

% IF NECESSARY, INTERCHANGE EQUATIONS TO MAKE THE LARGEST COEFFICIENT 
% BECOME THE PIVOTAL COEFFICIENT, 
if(imax ~= k) 

for j=k:jmax; % interchange coefficients 
atemp=a(k, j ); 
a(k,j)=a(imax,j); 
a(imax,j)=atemp; 
end; % interchange coefficients 
end; 

o, 

o 

% ELIMINATE X(K) FROM EQUATIONS (K+l) TO NEQN, FIRST TESTING FOR 
% EXCESSIVELY SMALL PIVOTAL COEFFICIENT (ASSOCIATED WITH A SINGULAR 
% OR VERY ILL-CONDITIONED MATRIX). 
if(abs(a(k,k)) < 1.Oe-5) 
iflag=l; 
return; 
end; 

for i=k+l:neqn; % each of remaining equations 
fact=a(i,k)/a(k,k); 

for j=k:jmax; % modify coefficients 
a(i,j)=a(i,j)-fact*a(k,j); 
end; % modify coefficient 
end; % each of remaining equations 

o 

o 

end; % eliminate each variable in turn 

o 

o 

% SOLVE THE EQUATIONS BY BACK SUBSTITUTION, FIRST TESTING 
% FOR AN EXCESSIVELY SMALL LAST DIAGONAL COEFFICIENT, 
if(abs(a(neqn,neqn)) < 1.0e-5) 
iflag=l; 
return; 
end; 

x(neqn)=a(neqn,jmax)/a(neqn,neqn); 

for i=neqn-l:-1:1; % then each unknown in turn backwards 
sum=a(i,jmax); 

for j=i+l:neqn; % sum products 
sum=sum-a(i,j)*x(j ) ; 
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end; % sum products 
x(i)=sum/a (i,i); 

end; % then each unknown in turn backwards 
return; 

end %function elimin 

The equivalent of the Fortran SHAREDDDATA2EQ is shareddata2eq, stored as the file shareddata2eq.m, 
is as follows. 


%%%- module shareddata2eq; 


% module STORING SHARED DATA. 


global xeend; if isempty(xeend), xeend=zeros(1,260); end; 
global yeend; if isempty(yeend), yeend=zeros(1,260); end; 
global xnode; if isempty(xnode), xnode=zeros(1,500); end; 
global ynode; if isempty(ynode), ynode=zeros(1,500); end; 
global xsend; if isempty(xsend ), xsend=zeros(1,250); end; 
global ysend; if isempty(ysend), ysend=zeros(1,250); end; 
global unmx; if isempty(unmx), unmx=zeros(250,3); end; 
global unmy; if isempty(unmy), unmy=zeros(250,3); end; 
global angstore; if isempty(angstore), angstore=zeros(1,260) 
global maxi; if isempty(maxi), maxl=0; end; 
global a; if isempty(a), a=zeros(1000,1001); end; 
global uv; if isempty(uv), uv=zeros(1,1000); end; 
global u; if isempty(u), u=zeros(1 , 500); end; 
global v; if isempty (v), v=zeros(1 , 500); end; 
global useg; if isempty(useg), useg=zeros(1 , 250); end; 
global vseg; if isempty(vseg), vseg=zeros(1,250); end; 
global arowx; if isempty(arowx ), arowx=zeros(250 , 6); end; 
global arowy; if isempty(arowy), arowy=zeros(250,6); end; 
global browx; if isempty(browx), browx=zeros(250 , 6); end; 
global browy; if isempty(browy), browy=zeros(250 , 6); end; 
global zg; if isempty(zg) A zg=zeros(1 , 8); end; 
global wg; if isempty(wg), wg=zeros(1 f 8); end; 
global egl; if isempty(egl), egl=zeros(1 f 8); end; 
global wgl; if isempty(wgl), wgl=zeros(1 f 8); end; 
global jacob; if isempty(jacob), jacob=0; end; 
global ungx; if isempty(ungx ), ungx=0; end; 


end; 
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global ungy; if isempty(ungy), ungy=0; end; 

global signnseg; if isempty(signnseg), signnseg~zeros(1,250); end; 
global sigsnseg; if isempty(sigsnseg), sigsnseg=zeros(1, 250); end; 


global 

uelem; if 

isempty(uelem), 

uelem=zeros(1,250); 

end; 

global 

velem; if 

isempty(velem), 

velem=zeros(1,250); 

end; 

global 

signn; if 

isempty(signn), 

signn=zeros (250,3); 

end; 

global 

sigss; if 

isempty(sigss), 

sigss=zeros (250,3); 

end; 

global 

sigsn; if 

isempty(sigsn), 

sigsn=zeros (250,3); 

end; 


global sige; if isempty(sige), sige=zeros(250,3); end; 
global tx; if isempty(tx), tx=zeros(1, 500); end; 
global ty; if isempty(ty), ty=zeros(1, 500); end; 
global tmx; if isempty(tmx), tmx=zeros(250,3); end; 
global tmy; if isempty(tmy), tmy=zeros(250,3); end; 
global fxseg; if isempty(fxseg), fxseg=zeros(1,250); end; 
global fyseg; if isempty(fyseg), fyseg=zeros(1,250); end; 
global sf; if isempty(sf), sf=zeros(3,8); end; 
global sd; if isempty(sd), sd=zeros(3,13); end; 
global sfl; if isempty (sfl) , sfl-zeros (4,3,8) ; end; 
global sdl; if isempty(sdl), sdl=zeros(4,3,8) ; end; 
global e; if isempty(e), e=0; end; 
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global nu; if isempty(nu), nu=0; end; 

global estore; if isempty(estore), estore=0; end; 

global akxx; if isempty(akxx), akxx=0; end; 

global akxy; if isempty(akxy), akxy=0; end; 

global akyx; if isempty(akyx), akyx=0; end; 

global akyy; if isempty(akyy), akyy=0; end; 

global bkxx; if isempty(bkxx), bkxx=0; end; 

global bkxy; if isempty(bkxy), bkxy=0; end; 

global bkyx; if isempty(bkyx), bkyx=0; end; 

global bkyy; if isempty(bkyy), bkyy=0; end; 

global bk2xx; if isempty(bk2xx), bk2xx=0; end; 

global bk2xy; if isempty(bk2xy), bk2xy=0; end; 

global bk2yx; if isempty(bk2yx), bk2yx=0; end; 

global bk2yy; if isempty(bk2yy), bk2yy=0; end; 

global nel; if isempty(nel), nel=0; end; 

global nnp; if isempty(nnp), nnp=0; end; 

global maxnel; if isempty(maxnel ), maxnel=0; end; 

global maxnnp; if isempty(maxnnp)^ maxnnp=0; end; 

global maxnb; if isempty(maxnb), maxnb=0; end; 

global neend; if isempty(neend), neend=0; end; 

global ngauss; if isempty(ngauss ), ngauss=0; end; 

global node; if isempty(node ), node=zeros(250 A 3); end; 

global ml; if isempty(ml), ml=zeros(1,250); end; 

global m3; if isempty(m3), m3=zeros(1,250); end; 

global nbound; if isempty(nbound), nbound=0; end; 

global nsegtot; if isempty(nsegtot), nsegtot=0; end; 

global nelb; if isempty(nelb), nelb=zeros(1,10); end; 

global nsegb; if isempty(nsegb), nsegb=zeros(1,10); end; 

global nbcu; if isempty(nbcu), nbcu=0; end; 

global nbcs; if isempty(nbcs), nbcs=0; end; 

global nbcm; if isempty(nbcm), nbcm=0; end; 

global nbct; if isempty(nbct), nbct=0; end; 

global ibce; if isempty(ibce), ibce=zeros(1,250); end; 

global ibcn; if isempty(ibcn), ibcn=zeros(1,500); end; 

global isegbc; if isempty(isegbc), isegbc=zeros(1,250); end; 

global isegend; if isempty(isegend), isegend=zeros(1,250); end; 

global isegelem; if isempty(isegelem), isegelem=zeros(1,250); end; 

global mfirst; if isempty (mfirst) , mfirst=zeros (1,250) ; end; 

global mlast; if isempty(mlast), mlast=zeros(1,250); end; 
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Appendix E: Matlab Version of Quadratic Boundary 
Element Program for Plane Elastic Problems 


global 

global 

global 

global 


nbcpc; 

maxnpc; 

nodepc; 

idirpc; 


if isempty(nbcpc)^ 
if isempty(maxnpc) 
if isempty(nodepc) 
if isempty(idirpc) 


nbcpc=0; end; 
f maxnpc=0; end; 

, nodepc=zeros(1 , 20); 
, idirpc=zeros(1 , 20); 


end; 

end; 
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Solutions to Problems - Part II 


Numerical answers are given, followed where appropriate by detailed notes on the solution procedure 
and comments on the results obtained. 


Chapter 5 

5.1 Indistinguishable when v = 0 (v* = v and E* = E ) . Of little practical significance because no 
known materials have Poissons ratios of zero, although some have values close to zero. 

5.2 An incompressible material has a Poissons ratio of v = 1/2. In plane strain v* = 1 and 
£* = 4E/3. The typical displacement and traction kernel functions become 

u ,«&>.</) = jM ,n (;) + 

Uxy(.P'l)=^r x t y 

T xx (p.q) = -^fA^ 

T xy <p l q) = -±t,t y % 


Unlike displacement-based finite element methods, incompressible plane strain presents no 
difficulties to the present boundary element formulation. 


5.3 


t x n x + t y n y = —k n {un x + vn y ) 

tyfi x - t x fiy = —k s (yn x - un y ) 


5.4 The body forces per unit volume in Equations 1.6 and 1.7 are X = 0 and Y = —pg , and the 
expressions for stresses satisfy these equations. Strains e xx and e xy are zero, from both the strain 
definitions of Equations 1.2 and 1.3, and in terms of stresses from Equations 1.20 and 1.23. From 
Equations 1.2 


_ pay 

e yy ~ Z 7 * (1 — v ) 


and from Equation 1.21 


i r * \ pay sm *2 > v 

e yy =7^0yy °xx) =-^r(l-v 2 ) 


The expressions for e yy are identical, and all the relevant equations are satisfied. 


188 


Download free eBooks at bookboon.com 




Boundary Element Methods for Engineers: 
Part II: Plane Elastic Problems 


Solutions to Problems - Part 2 


5.5 The body forces per unit volume in Equations 1.6 and 1.7 are X = pa> 2 x and Y = 0, and the 
expressions for stresses satisfy these equations. Strains e yy and e xy are zero, from both the strain 
definitions of Equations 1.2 and 1.3, and in terms of stresses from Equations 1.21 and 1.23. From 
Equations 1.2 


e 


XX 


2 2 

porx z 

2 £* 


(1 - V* 2 ) 


and from Equation 1.20 


1 , . pa) 2 x 2 , _ 

= ^Oyy -v'ir a ) =-) 


The expressions for are identical, and all the relevant equations are satisfied. 


Chapter 6 

6.1 Answer: maximum hoop stress concentration factor 1.005. The analytical solution for a hole 
in an infinite plate (/f -> oo) is 1 (Equations 6.8). For a plate width 20 times the hole diameter 
(/f = 20), which ignores the square corners of the plate, it is 1.005, identical to four significant 
figures to that computed. But, with compressive radial stress (the pressure) at the surface of the 
hole, the hoop stress is not the largest stress to be considered there. The maximum von Mises 
equivalent stress concentration factor (relative to the applied pressure) is 1.74. Results were 
obtained using 4 elements per side of the outer boundary, 6 elements per quadrant of the hole. 

6.2 Answers: at the inner surface, hoop stress 0.231, radial displacement 0.00240. At the outer 
surface, radial stress -0.539, hoop stress -0.231. Stress values are averaged over element end 
nodes and midside nodes: for example, typical hoop stress values at the inner surface are 0.2313 
at end nodes, 0.2309 at midside nodes (a difference which is not significantly reduced by mesh 
refinement). Results were obtained using 4 elements per quadrant on both boundaries. 

The general analytical solutions for stresses in a thick-walled cylinder are given by Equations 

6.5. In this problem the boundary conditions are o rr — — p at r — r x and e ee — 0 at 
r = r 2 > from which 

_ —p(l+v*)r! 2 _ p(l—v*)r 2 r| 

(1—v*)r|+(l+v*)r 2 5 (1—v*)r| +(l+v*)r 2 

At the inner surfaces = ?y and the hoop stress is 
= p[(l-v*)fr 2 -(l+v*)] 

(i— v *)ff2 +( i +v .) 
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where, K = ^ = 2,p = 1 and v* = - = 0.42857, giving a ee = 0.2308. 

At the outer surface, r — r 2 and the radial and hoop stresses are 

a rr = (1 _ v . ) ^2+ (1+v » ) = -0.5384, Ggg = V*G rr = -0.2308 

Radial displacement at the inner surface can be found from the hoop strain with the aid of 
Equation 6.9 

^7 (^00 V ®rr ) 

Where a ee — 0.2308, cr r = — 1, and E* = _^ 2 = 1098.9, giving u r — 0.002402. 

6.3 Answer: 3.12. Result obtained using 4 elements per edge of the plate on edges without the notch, 
4 elements per quadrant in the notch, and 10 elements on each of the segments adjacent to the 
notch. No significant change in the result was found when element sizes on these two segments 
were graded to give smaller elements near the notch. 


Answer: 7.81. Result obtained using 4 elements per outer edge of the plate, 20 per straight side 
of the square hole, and 8 per quadrant of the rounded corners. 
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6.5 Answer: 0.681, a friction coefficient of more than 0.5 would be required to prevent slipping. 
The friction requirement is investigated by finding the ratio between shear and normal stresses 
(SIGSN and SIGNN) at nodes on the base of the pedestal, giving the minimum coefficient to 
prevent slipping. In this case the value of 0.5 is exceeded at points more than about 0.07 m 
from the central line of symmetry (less than about 0.08 m from the edges of the pedestal base). 
Results were obtained using 4 elements on each of the top and sloping sides of the pedestal, 
and 8 on the base. The boundary conditions are zero displacements on the base, and a normal 
stress of -1 on the top of the pedestal. Such contact problems, involving slippage with friction 
and possible local loss of contact, can be solved by boundary element methods, but are non¬ 
linear in the sense that the boundary conditions change with the loading, and require more 
sophisticated programming. 

6.6 Answers: stiffness 142 MN/m, maximum stress (SIGSS) 190 MN/m 2 at the fillets, stress (SIGSS) 
variation over the central 30 mm of length is from 89.3 to 89.7 MN/m 2 , compared with the 
nominal stress of 90 MN/m 2 . A state of plane stress maybe assumed. Results were obtained using 
4 elements per straight edge at the ends of the test piece, 4 elements on each of the fillets, and 10 
elements along each side of the central section (to provide more detail there, and to have element 
nodes at the ends of the central 30 mm). Point constraints were applied to the corners at the left 
hand end of the test piece. The stiffness is estimated from the computed relative displacement 
between the nodes at the centres of the two ends: 0.380 X 10 -4 m. The tensile force in the test 
piece is 60 X 10 6 X 0.03 X 0.003 = 5400 N (stress times width times thickness), giving a 
stiffness of 5400/0.380 X 10 -4 N/m or 142 MN/m. 

6.7 To apply pure shear, balancing shear stresses must be applied to all four edges of the plate: 
for example, <7 sn = +1 on BC and DA, o sn = —1 on CD and AB (Figure 6.7). The maximum 
computed hoop stresses at the hole are 4.03, occurring at four angular positions at 45° to the 
x and y axes. The analytical solution for a hole at the centre of an infinite plate in pure shear 
gives a value of 4 times the applied shear stress. Results were obtained using 4 elements per side 
of the outer boundary, 6 elements per quadrant of the hole. 

6.8 Answer: 5.03 (hoop stress concentration factor at both sides of the hole). Rather than create 
the mesh data for an elliptical hole boundary, the same geometric data as in Problems 6.1 and 6.7 can 
used, but with they coordinates of the corners of the outer boundary doubled in magnitude. Then in 
subprogram MESHQ all YNODE values can be scaled by a factor of 0.5 immediately before finding 
the maximum dimension of the domain. This creates the elliptical shape of hole and brings the outer 
boundary shape back to square. Results were obtained using 4 elements per side of the outer boundary, 
6 elements per quadrant of the hole. The exact solution for an elliptical hole with a ratio of major to 
minor axes of 2, under tension in an infinite plate, is 5. 
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6.9 Answers: (a) 3.11 at the outer sides of the holes (y = 0,x = ±0.075), 3.06 at the inner sides 
(y = 0,x = ±0.025). (b) 2.74 (at x = ±0.05 ,y = ±0.025). Results were obtained using 4 
elements per side of the outer boundary, 6 elements per quadrant of each of the holes. An 
example of a solution domain with three boundaries. 

6.10 Answer: 745 MN/m 2 at the bottom of the groove. The position of the maximum equivalent stress 
is as expected. With a radius ratio of K = 2.5, the hoop stresses at the inner and outer surface 
of the cylinder in the absence of the groove are 

o gg = p = 166 MN/m 2 MN/m 2 and age = = 45.7 MN/m 2 

These figures are useful for checking the computed results: hoop stresses at the inner surface 
remote from the groove should be reasonably close to these (within about 1-2%, say). The 
choice of mesh in this problem requires some care. Experience shows that 4 elements per 
quadrant should be sufficient on the outer boundary, and perhaps 6 elements per quadrant 
in the semi-circular groove. On the inner cylindrical surface it is important to have elements 
near the groove which are similar in length to the very small elements in the groove. Using 
elements of uniform size to satisfy this requirement leads to an unnecessarily large number of 
elements. A reasonable compromise was found to be 16 elements on each of two nearly semi¬ 
circular segments (that is, about 8 per quadrant), with an element length ratio (as explained in 
Sections 3.1.3 and 4.1.3) of. In other words, the constant ratio between the lengths of successive 
elements moving away from the groove was 1.3. The total number of elements used was 60. In 
addition to checking inner surface hoop stresses remote from the groove against the analytical 
values above, another useful test is for the tangential stresses (SIGSS) in the elements adjacent 
to the corner between the groove and the inside surface of the cylinder: at the common node 
they should both be close to the applied pressure. The computed hoop stress at the bottom of 
the groove is 678 MN/m 2 , 4.08 times the hoop stress at the inner surface in the absence of the 
groove. Combined with a radial stress of -120 MN/m 2 , this gives a von Mises equivalent stress 
of 745 MN/m 2 , some 6.2 times the applied pressure in magnitude, also some 3.00 times the 
maximum equivalent stress in the absence of the groove. 
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6.11 Answer: 5.23. The problem does not state whether the fault in manufacture caused the cylinder 
to distort after machining, in which case the wall thickness would be constant around the 
circumference, or both inner and outer surface were made slightly elliptical. The effect on the 
result is unlikely to be significant at such a low level of eccentricity, and for convenience the 
latter is assumed. The change in geometry is effected by modifying the original mesh to scale all 
the YNODE values by a factor of 0.99 in subprogram MESHQ immediately before finding the 
maximum dimension of the domain. Results were obtained using 4 elements per quadrant on both 
boundaries. The result for maximum hoop stress in the circular cylinder is 5 times the internal 
pressure. So, a 1% eccentricity of the cylinder gives a 4.6% increase in the maximum hoop stress. 

6.12 Answer: 299 Nm. The problem did not specify whether plane stress or plane strain is to be 
assumed. Because the state of deformation in the rubber bush is one of shearing, plane stress 
and plane strain give the same results. The problem calls for torque to be computed (via shear 
stress on the inner boundary) for a given rotation of this boundary. This rotation displacement 
boundary condition cannot be applied using the present program. On the other hand, a uniform 
shear stress can be applied to the inner boundary: a value of 1 N/m 2 was used. Results were 
obtained using 4 elements per quadrant on both the inner and outer circular boundaries. The 
computed value of tangential displacement at the inner boundary was 0.15750 X 10 -9 m. At 
a radius of 0.01 m, this corresponds to an angular displacement of 0.15750 X 10 -7 radians. 
The shear stress required to produce the required rotation of 0.1 radians is therefore 

shear stress = 1 X - ——^ 7 = 6.3492 X 10 6 N/m 2 

0.15750 XlO -7 

This acts over a surface area of radius 0.01 m and length 0.075 m at a distance of 0.01 m from 
the axis of the shaft. Hence, the torque is 

torque = 6.3492 X 10 6 X 2 X n X 0.01 2 X 0.075 = 299 Nm 
The problem can also be solved analytically, giving a result of 299.2 Nm. 
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